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ABSTRACT 


The  work  described  is  a  continuation  of  that  reported  previously 
in  AEDC-TN-6l”0^  and  AEDC-TDR-62-131.  In  particular,  the  present 

report  contains  the  basis  far  She  ■athad  dem.!  lbl!a"n!l . ABWTPn  60  tSl  - 

for  the  precise  numerical  calculation  of  one -dimensional  nonequilibrium 

flow  of  a  complex  gas  mixture  through  a  nozzle K  The  contributions -ef - 

tha  two  above  reports  are  »'eviewed~'bi  ief  ly  in- the  4f»t«>«»d«ct.ion. 


A  transformation  is  shown  to  exist  when  the  flow  is  close  to  equi¬ 
librium  that  transforms  the  chemical  rate  equations  into  a  simpler 
system  of  uncoupled  rate  equations.  Each  of  the  transformed  chemical 
rate  equations  is  similar  in  form  to  a  vibrational  rate  equation.  The 
coefficients  appearing  in  the  transformed  equations  are  given  by  solu¬ 
tions  of  an  eigenvalue  problem.  These  coefficients  are  shown  to  be 
real  and  positive.  JAn  example  is  given  that  illustrates  the  transfer- 
matlon  for  a  simple  gas  model  and  shows  that  the  transformation  behaves 
correctly  in  the  frozen  and  equilibrium  limits. 'y 

Because  of  their  simple  mathematical  structure,  the  vibrational 
and  transformed  chemical  rate  equations  can  be  subjected  to  numerical 
analysis.  Thus,  a  detailed  analysis  of  both  Runge-Kutta  and  predictor- 
corrector  integration  procedures  is  performed  in  which  truncation  error, 
stability,  and  computation  speed  are  examined.  The  analysis  of  Runge- 
Kutta  procedures  demonstrates  that  the  greater  the  order  of  the  proce¬ 
dure,  the  more  stable  it  is.  High-order  Runge-Kutta  procedures,  unfer-~ 
_tunafee-iy>  also  have  reduced  computation  speed.  On  an  overall  basis, 
the  commonly  employed  fourth-order  procedure  still  appears  to  be  the 
most  suitable  Runge-Kutta  procedure  for  the  integration  of  the  chemical 
and  vibrational  rate  equations.  The  basis  for  bhK improved  technique 

for  controlling  the  Runge-Kutta  integration  step  size,  girew  In - 

AEDC-TDR-69-131,  is  also  established. 

Vfhen  Adame  predictor-corrector  procedures  are  similarly  analyzed, 
the  greater  the  order  of  the  procedure  the  less  stable  it  is.  Further¬ 
more,  the  overall  Adams  procedure  is  considerably  less  stable  than  is 
the  Adams  corrector  formula.  Tlie  common  assumption  that  stability  is 
determined  primarily  by  the  corrector  is  thus  not  valid.  A  fourth- 
order  predictor-corrector  procedure  is  also  given  that  is  stable  for 
a  larger  integration  step  size  than  the  fourth-order  Adams  procedure. 
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stability  parameter  for  Runge-Kutta  Integration 
procedures 

coefficients  In  the  transformation  (l9a) 

stability  parameter  for  predictor  formula 

stability  parameter  for  corrector  formula 

stability  parameter  for  overall  predictor-corrector 
procedure 

characteristic  polynomial  for  predictor  formula 

characteristic  polynomial  for  corrector  formula 

characteristic  polynomial  for  overall  predictor - 
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local  equilibrium  vibrational  energy  corresponding  to 
static  temperature  T 
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defined  by  equation  (l-l6) 
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integration  step  size 

defined  by  equations  (1-13)  and  (l-15a) 

equilibrium  constant  for  reaction  r 

see  predictor-corrector  formulas  (48)  and  (4-9) 

backward  rate  constant  for  reaction  r 


forward  rate  constant  for  reaction  r 
see  equation  (6b) 

total  number  of  reacting  chemical  species 


number  of  chemical  reactions 


number  of  chemical  species  minus  total  number  of 
components  and  polyatomic  inert  species 

number  of  chemical  species  plus  number  of  catalytic 
bodies  that  consist  of  a  linear  combination  of 
other  species 

mole-mass  ratio  of  species  i  ;  that  is,  number  of 
moles  of  species  i  per  unit  mass  of  fluid 

local  equilibrium  value  of  n^  based  on  the  local 
static  density  and  temperature 

defined  by  equation  (l-8a) 


order  of  the  truncation  procedure 

generalized  mole-mass  ratio  (see  equation  (l9a)) 


static  temperature 
flow  speed 


distance  along  nozzle  axis 
variables  in  equation  (2a) 
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*(x^4ll) 
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ri 


rl 


k/ 


estinate  for  y(x^+h)  given  by  predictor  foraula  (48) 
estimate  for  y'(x^+h)  used  in  corrector  formula  (49) 
constants  in  the  Range -Kutta  integration  procedure 
see  equation  (7) 

forward  stoichiometric  coefficient  for  species  i  in 
reaction  r 

backweurd  stoichiometric  coefficient  for  species  i  in 
reaction  r 

see  equation  (4la) 

Kronecker  delta 

small  constant  (see  equation  (2b)) 


6 

r 


defined  by  equation  (6c) 


K 


"U 


unknovm  in  the  characteristic  polynomial  given  by 
equation  (24b) 

roots  of  the  characteristic  polynomial  given  by 
equation  (24b) 

defined  by  equation  (l6) 
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defined  by  equation  (l8) 


rl 


stoichiometric  coefficients  (see  equation  (6a)) 


P 


mass  density  of  fluid 


P 

Pi 


unknown  in  the  stability  characteristic  polynomial 
roots  of  the  stability  characteristic  polynomial 
vibrational  relaxation  time  for  species  i 


AeOC.TDR  <s.ia 


SUKRSCMm  AND  SUMCRim 

(  )*  equilibrium  reference  value 

(  ) '  perturbation  value 

(  ) '  denotes  differentiation  with  respect  to  x 

e  value  at  local  equilibrium  condition 

1  value  for  species  1 

indices 

In  initial  ralue 

o  value  at  the  beginning  of  an  integration  step 

r  value  for  reaction  r 
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1.0  INTRODUCTION 

This  is  the  third  In  a  series  of  reports  concerning  the  calculation 
of  the  reacting  flow  of  a  cooplex  gas  in  the  nozzle  of  a  hypersonic 
wind  tunnel  at  high  stagnation  enthalpy.  The  first  of  these  reports 
(Vincent!  (196I),  hereafter  referred  to  as  Part  l)  describes  a  five- 
species  model  for  air,  governed  by  eight  chemical  kinetic  reactions. 

A  method  is  described,  in  particulcur,  for  the  numerical  calculation  of 
the  one -dimensional  nonequilibrium  flow  of  this  gas  through  a  hypersonic 
nozzle.  A  specific  calculation  ccurrled  out  on  an  IBM  709  con^ter,  how¬ 
ever,  revealed  that  the  method  required  too  much  conqniter  time  to  be 
practical  for  engineering  purposes. 

The  second  report  (Emanuel  and  Vlncentl  (1962),  hereafter  referred 
to  as  Part  II)  describes  an  Improved  method  for  calculating  nonequl- 
llbrlum  nozzle  flews.  This  method  Is  relatively  slnqple,  niimerlcally 
accurate,  and  does  not  require  an  excessive  amount  of  machine  time.  It 
can  be  summarized  as  follows: 

(a)  An  equilibrium  solution  for  the  nozzle  flow  is  first  obtained. 
To  facilitate  this,  a  new,  simple  procedure  for  computing  an 
Inexact  but  numerlcsilly  accurate  equilibrium  solution  was 
given. 

(b)  Equilibrium  Initial  conditions  are  used  to  start  the  nonequl- 
llbrlvun  calculation.  This  simple  method  for  starting  the 
nonequillbriiun  calculation  replaces  the  conpllcated  pertur¬ 
bation  schemes  frequently  employed. 

(c)  Forward  Integration  of  the  nonequlllbrlvm  equations  then  pro¬ 
ceeds  from  the  equilibrium  Initial  point.  This  point  Is 
chosen  somewhat  upstream  of  the  nozzle  location  at  which  the 
chemistry  first  departs  appreciably  from  the  equilibrium- 
flow  solution. 

With  this  method,  the  numerical  solution  Is  Insensitive  to  the  location 
of  the  equilibrium  Initial  point,  provided  that  It  is  chosen  as  described 
In  (c)  above. 

The  above  method  requires  that  the  nonequilibrium  calculation 
Includes  a  portion  of  the  nozzle  in  which  the  flow  is  close  to  equlllbrliun. 
Although  the  extent  of  this  region  can  be  minimized,  a  very  small  step 
size  Is  nevertheless  necessary  for  the  Integration  therein.  If  too 
large  a  step  size  is  used  in  this  near-equlllbrlum  region,  then  the 
Integration  procedure  becomes  unstable  and  the  numerical  solution  di¬ 
verges  from  the  correct  solution.  Thus,  to  avoid  an  excessive  amoxmt 
of  conputer  time,  the  Integration  procedure  must  allow  the  step  size  to 
vary  In  a  manner  that  maximizes  Its  size  but  still  keeps  the  procedure 
stable.  Part  II  gave,  for  a  Runge-Kutta  Integration  procedure,  a  simple 
method  of  controlling  step  size  that  approximately  satisfies  the  fore¬ 
going  criterion. 
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The  present  report  contains  the  analytical  basis  for  the  numerical 
method  described  in  Psort  II  and  outlined  above.  Therefore,  some 
acquaintance  by  the  reader  with  the  system  of  equations  of  Part  II  Is 
desirable.  We  shall  first  shew  that  a  system  of  chemical  and  vibra¬ 
tional  rate  equations  Is  equivalent,  for  purposes  of  numerical  analysis, 
to  a  slniple,  so-called  "stiff"  equation.  The  remainder  of  the  report 
Is  then  concerned  with  the  numerical  analysis  of  this  equation. 

The  author  wishes  to  express  his  sincere  gratitude  to  Professor 
Valter  G.  Vlncentl  for  his  guidemce  In  this  work.  Thanks  are  also  due 
to  Mrs.  Llta  Emanuel  for  editorial  assistance. 


2.0  PRELIMiNARY  REMARKS 


As  Curtiss  and  Hlrschf elder  (1952)  have  pointed  out,  a  type  of 
differential  equation  occurs  In  many  physical  problems  that  is  extremely 
difficult  to  solve  numerically.  They  designate  this  type  of  equation 
as  "stiff"  since  a  typical  example  Is  the  equation  that  describes  the 
motion  of  a  simple  mechanical  system  with  a  stiff  spring.  We  shall  here 
define  a  stiff  equation  as  any  ordinary  or  peurtial  differential  equation 
in  which  the  hlghest-order  derivative  Is  multiplied  by  a  small  i>arameter. 
Problems  governed  by  this  type  of  equation  are  also  referred  to  as 
singular  perturbation  problems.  This  designation,  however,  is  correct 
only  in  the  limit  as  the  small  parameter  goes  to  zero.  Since  our  inter¬ 
est  here  Is  in  the  numerical  analysis  of  this  type  of  equation,  and  not 
in  the  limiting  process,  the  terminology  of  Curtiss  and  Hlrschf elder 
will  be  used. 

In  a  near-equlllbrlum  flow,  each  chemical  and  vibration  rate  equa¬ 
tion  Involves  a  small  parameter  that  multiplies  the  derivative.  Thus, 
in  a  near-equlllbrlum  flow,  each  rate  equation  is  stiff.  For  example, 
uTj^  is  the  small  parameter  in  the  vibrational  rate  equations 


(1) 


which  applies  to  steady  one -dimensional  flow.  Here  e^  euid  e^ 

are,  respectively,  the  actual  vibrational  energy  and  the  local  equiflb- 
rltun  vibrational  energy  per  mole  of  species  1  .  The  flow  speed  Is 
denoted  by  u  ,  vibrational  relaxation  time  of  species  1  by  t.  ,  and 
x  is  the  distance  along  the  nozzle  axis.  (For  additional  details  see 
either  Parts  I  or  II.) 


Our  object  Is  the  ansQysls  of  stiff  rate  equations,  both  chemical 
and  vibrational.  Because  of  their  complexity,  however,  the  chemical 
rate  equations  will  not  be  studied  directly.  Instead,  they  will  first 
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be  transformed  into  a  system  of  equations  of  the  same  form  as  the  vibra¬ 
tional  rate  equations  (l).  It  will  therefore  not  be  necessary  to  dis¬ 
tinguish  between  the  two  types  of  equations  in  the  analysis. 

To  derive  the  transformation,  which  is  given  in  detail  in  Chapter 
3.0  and  illustrated  by  an  example  in  Appendix  II,  the  chemical  rate 
equations  are  first  linearized  about  an  equilibrium  reference  state. 

As  a  result,  the  transformation  is  valid  only  for  a  near-equlllbrlum 
flow.  In  general,  the  reference  state  is  variable,  as  in  nozzi^e  flow, 
and  the  transformation  is  then  only  a  local  one.  In  other  words,  the 
transformation  is  vsilld  only  over  a  limited  region  of  the  flow.  A 
sequence  of  transformations  may  thus  be  required  to  represent  the  entire 
near-equlllbrlum  region.  Our  purpose  in  deriving  the  transformation, 
however,  is  to  be  able  to  analyze  numerically  the  behavior  of  the  rate 
equations  when  the  flow  is  close  to  equilibrium.  Since  this  analysis 
is  also  of  a  local  nature,  the  limited  region  of  validity  of  the  trans¬ 
formation  is  of  no  conseqiwnce. 

Inasmuch  as  the  difficulties  of  integration  of  the  nonequilibrium 
equations  arise  solely  out  of  the  stiffness  of  the  rate  equation, 
certain  assumptions  may  be  introduced  to  simplify  the  analysis.  These 
assuiiq>tions,  however,  in  no  way  alter  or  reduce  the  numerical  diffi¬ 
culties  from  those  encountered  when  the  original  rate  equations  are 
Integrated.  To  better  understand  these  assumptions,  we  first  note  that 
the  vibrational  rate  equations  (l)  are  coupled  only  through  the  depend¬ 
ence  of  e„  and  ut  on  the  density  and  temperature  of  the  flow. 

The  same  is  also  true  of  the  transformed  chemical  rate  equations  (see 
equations  (20)),  which  for  steady  one-dimensional  flow  have  the  same 
structure  as  the  vibrational  equations.  Thus,  there  will  be  no  loss 
of  generality  if  the  following  is  assxaned  (eifter  the  chemical  rate 
equations  have  been  transformed  into  their  simplified  form); 

(a)  All  of  the  rate  equations  are  uncoupled.  This  is  equivalent 
to  assuming  that  the  density,  temperature  and  flow  speed  are 
known  functions  of  the  indei)endent  variable  x  ,  since  any 
coupling  in  the  rate  equations  occurs  through  these  variables. 
More  specifically,  when  the  equilibrium-flow  solution  is 
known,  these  variables  may  be  approximated  by  their  equilib¬ 
rium-flow  values.  This  assumption  in  no  way  alters  the  struc¬ 
ture  of  the  rate  equations,  and  only  slightly  modifies  the 
magnitudes  of  the  small  parameters.  Hence,  as  stated  earlier, 
no  loss  of  generality  ensues  from  this  decoupling  assumption. 

(b)  In  place  of  a  system  of  uncoupled  rate  equations,  it  is 
sufficient  to  consider  only  one  such  equation.  This  equation 
is  taken  to  be  that  containing  the  smallest  parameter.  When 
this  assumption  is  examined  a  posteriori,  it  is  readily  found 
to  be  Justified. 
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An  Inportant  conseq.uence  of  the  t;ranafomation  referred  to  earlier 
is  that  the  enall  paraoetera  in  the  transforeed  ehenieal  rate  equationa 
are  all  poaitive.  (Appendix  I  eontaina  the  proof  for  thia  atatenent.) 
Thia  fact 4  in  conjunction  with  the  foregoing  aaauiiqptiona,  inpliea  that 
ve  need  exaaiine  only  one  equation  with  the  aauBe  form  as  the  vibrational 
rate  equationa  (l).  Under  these  assumptions,  the  stiff  equations  nay 
be  represented  by 


dx 


(2a) 


and  euiy  numerical  analysis  that  is  valid  for  this  equation  will  be 
equally  valid  for  the  original  system  of  chemical  and  vibrational  rate 
equations  of  Farts  I  and  II.  In  equation  (2a),  e(x)  Is  a  small 
positive  function,  while  0(x)  Is  an  curbltrary  function.  Equation  (2a) 
can  be  further  simplified,  however,  still  without  any  loss  of  general¬ 
ity  Insofar  as  the  numerical  analysis  Is  concerned.  To  this  end,  we 
now  also  assume  the  following: 

(c)  The  functions  e(x)  and  G(x)  are  given  by 

€(x)  =  small  positive  constant,  ^ 

G(x)  *  X  .  J 

These  simplifications  sure  Justified  by  first  expanding  e 
and  G  In  Taylor  series  about  any  point  x'  .  The  variables 
X,  y  are  then  transformed  so  that  the  numerator  of  the 
rl^t-hand  side  of  equation  (2a)  has  (in  terms  of  the  new 
variables)  the  form  y+x  +(con8tant)  x^+  ...  .  Finally, 

only  the  lowest  order  terms  are  retained  in  the  right-hand 
side  of  equation  (2a),  which  leads  to  equations  (2b).  This 
last  step  Is  consistent  with  the  local  nature  of  the  initial 
transformation  of  the  chemical  equations  as  well  as  with  the 
subsequent  numerical  analysis. 

Because  of  their  simplicity  as  compared  with  the  original  rate 
equations,  equations  (2)  are  readily  analyzed.  From  this  suialysis,  two 
Important  facts  become  apparent.  First,  at  the  outset  of  the  nonequl- 
llbrlum  calculation  a  small  Integration  step  size  is  necessary.  This 
requirement  Is  a  direct  consequence  of  the  magnitude  of  the  truncation 
error  when  the  Integration  starts  from  an  equilibrium  condition.  Second, 
the  Integration  step  size  can  Increase  appreciably  after  a  short  ficti¬ 
tious  transient  region  Is  completed  wherein  the  derivatives  of  the  vari¬ 
ables  change  from  their  frozen  values  to  approximately  their  equilibrium 
values.  After  this  transient  region,  the  size  that  the  Integration  step 
can  attain  depends  on  the  specific  Integration  procedure.  For  this 
reason,  the  two  most  commonly  used  procedures  for  numerical  Integration, 
Runge-Kutta  (Chapter  5*0)  and  predictor -corrector  (Chapter  6.0),  are 
analyzed.  This  analysis,  it  should  be  noted.  Is  applicable  only  to  the 
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Integration  of  the  nonequilibrium  equations  when  the  flow  is  close  to 
equilibrium.  Once  the  integration  variables  have  begun  to  diverge  from 
their  local  equilibrium  values,  the  nonequilibrium  equations  are  easily 
Integrated,  and  any  standard  procedure  is  satisfactory.  The  conqparatlve 
evaluation  of  the  two  procedures,  given  in  Chapter  7*0,  reveals  that  a 
particular  fourth-order  predlctor-coirector  procedure  appears  to  be  the 
most  suitable  for  the  integration.  It  is  not,  however,  appreciably 
superior  to  a  fourth-order  Runge-Kutta  procedure.  Because  this  latter 
procedure  is  considerably  sloqpler  to  code  for  a  computer  than  are  the 
other  procedures,  it  was  utilized  for  all  the  numerical  integrations 
given  in  Part  II. 

The  author,  in  an  extensive  search  of  the  literature,  could  locate 
only  two  published  articles  dealing  directly  with  the  numerical  inte¬ 
gration  of  stiff  equations.  The  first  is  the  previously  mentioned 
paper  by  Curtiss  and  Hirschfelder  (1952),  which  recommends  the  use  of 
a  certain,  simple  Integration  procedure.  Inasmuch  as  their  analysis 
is  misleading,  this  procedure  is  discussed  in  Section  6.3. 

The  second  article  is  by  Certaine  and  is  contained  in  the  volume 
edited  by  Ralston  and  Wilf  (1960).  Hamming  (1962)  also  presents 
Certaine 's  procedure  in  a  simplified  form.  The  procedure  is  as  follows: 

(a)  First,  add  y/c  to  both  sides  of  equation  (2a)  and  then 
multiply  by  exp(x/6)dx  ,  thereby  obtaining 

exp(x/€)dy  +  y  exp(x/€)  ^  =  -  7  exp(x/€)G(x)dx  .  (3a) 


(b)  If  e  is  assumed  constant  over  an  Interval,  say  from  x^ 
to  Xg  ,  then  equation  (3a)  may  be  written  as 

d^y  exp(x/€)^  “  “  c  exp(x/€)G(x)dx  .  (3b) 

(c)  Finally,  Integrate  equation  (3b)  from  x^  to  Xg  and  then 
rearrange  the  results,  thereby  obtaining 

yg  =  -7  exp(-Xg/e)  f  exp(x/€)G(x)dx  ,  (3c) 

'  *1 

where  the  quadrature  on  the  right  is  readily  obtained  by 
standard  numerical  techniques. 
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The  principal  objection  to  Certalne's  procedure  Is  that  It  cannot 
be  used  to  Integrate  the  chenlcal  rate  equations  given  in  Parts  I  or  II 
unless  their  right-hand  sides  are  first  linearized  and  the  snsdl  para¬ 
meters  explicitly  exhibited.  A  second  difficulty  Inherent  In  the  pro¬ 
cedure  Is  the  choice  of  the  appropriate  Integration  step  size  X2  -  . 

In  partlculeur,  the  step  size  must  be  sufficiently  small  that  the  small 
parameters  core  nearly  constant  over  the  Integration  Interval.  Thus, 
the  rate  of  change  of  magnitude  of  each  small  parameter  must  be  period¬ 
ically  checked  to  see  whether  the  step  size  should  be  altered. 

3.0  CANONICAL  FORM  FOR  THE  CHEMICAL  RATE  EQUATIONS 

In  this  chapter  we  shall  show  that  when  the  flow  Is  close  to  equi¬ 
librium,  the  chemical  rate  equations  can  be  transformed  Into  a  system 
of  stiff  equations  having  the  scune  form  as  the  vibrational  rate 
equations.  While  this  transformation  Is  not  directly  useful  for  solv¬ 
ing  problems.  It  does  Justify  the  use  of  equations  (2),  which  are  basic 
to  the  subsequent  analysis.  The  particular  transformation  to  be  devel¬ 
oped  is  not  limited  to  steady  one -dimensional  flow,  but  Is  valid  for 
unsteady  three-dimensional  flow  as  well. 

The  situation  here  is  similar  to  that  of  a  conservative  mechanical 
system  that  performs  small  oscillations  about  a  fixed  equilibrium 
configuration.  Many  textbooks,  such  as  Goldstein  (19^),  show  that  a 
linear  transformation  can  be  introduced  to  transform  the  equations  of 
motion  Into  a  form  referred  to  as  canonical  because  of  Its  simple 
mathematical  structure.  Since  this  procedure  will  be  followed  approx¬ 
imately  in  our  derivation  of  the  canonical  form  for  the  chemical  rate 
equations,  the  general  procedure  for  a  vibrating  mechanical  system  will 
first  be  summarized  as  follows: 

(a)  The  potential  and  kinetic  energies  are  each  expanded  In  a 
Taylor  series  about  the  equilibrium  configuration.  Only  the 
first  approximation  to  the  potential  euid  kinetic  energies  Is 
retained.  Consequently,  the  equations  of  motion  are  linear, 
second -order  differential  equations. 

(b)  A  principal -axis  (or  normal -coordinate)  transformation  is 
Introduced  that  diagonalizes  the  coefficients  of  the  poten¬ 
tial  and  kinetic  energies  and  thereby  transforms  the  equa¬ 
tions  of  motion  Into  their  canonical  form. 

(c)  An  eigenvalue  problem  Is  solved  for  the  elgenfrequencles  in 
order  to  find  the  coefficients  In  the  principal -axis 
transformation.  Because  the  first  approximation  to  the 
potential  energy  is  positive  definite,  all  the  elgenfre¬ 
quencles  can  be  shown  to  be  real  and  positive. 
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In  parallel  with  the  foregoing,  the  right-hand  side  of  the  chemical 
rate  equations  will  first  be  linearized  about  an  equilibrium  reference 
state  in  a  manner  similar  to  that  given  by  Vincent!  (1959)-  The  refer¬ 
ence  state  for  a  slightly  varying  flew  may  be  chosen  as  the  uniform 
undlstxirbed  (or  free-stream)  condition.  When  the  flow  has  large  vari¬ 
ations,  as  in  a  nozzle,  the  reference  state  Is  taken  as  the  nonconstsuit 
equilibrium-flow  solution.  In  general,  we  shall  consider  the  reference 
state  to  be  a  function  of  position  and  time.  A  linear  transformation 
of  the  quantities  that  specify  the  chemical  composition  Is  then 
Introduced.  This  transformation,  when  subjected  to  the  condition  that 
the  rate  equations  trsuisform  Into  their  canonical  form,  results  In  an 
eigenvalue  problem.  Each  eigenvalue  Is  the  reciprocal  of  the  small 
parameter  In  a  canonical  rate  equation.  Because  of  Its  complexity  and 
length,  the  proof  that  all  eigenvalues  are  real  and  positive  Is  given 
In  Appendix  I. 

3.1  CHEMICAL  RATE  EQUATIONS 

The  following  four  quantities,  necessary  In  the  subsequent  analy¬ 
sis,  are  here  defined  as  follows  (cf.  Part  II): 

=  total  number  of  reacting  chemical  species, 

Ng  =  number  of  chemical  reactions, 

Nj  =  number  of  chemical  species  minus  total  number  of  components 
and  polyatomic  Inert  species, 

Nj^  =  number  of  chemical  species  plus  number  of  catalytic  bodies 
that  consist  of  a  linear  combination  of  other  species. 

As  an  Illustration,  given  also  In  Part  II,  consider  a  model  of  air  con¬ 
sisting  of  0,  N,  NO,  I^,  O2  ,  and  A  and  having  eight  chemical  reac¬ 
tions  to  specify  the  chemistry.  The  number  of  species  N-j^  is  then 
six,  and  the  number  of  reactions  Is  eigdit.  Since  there  are  three 
atomic  species  (one  of  which  is  Inert),  N3  is  three.  The  quantity 
N4  Is  Introduced  In  order  to  simplify  the  treatment  of  catalytic 
bodies  made  up  of  a  group  of  species,  such  as  N  +  Ng  +  NO  .  Thus,  a 
fictitious  species  M  ,  which  represents  any  atom  or  molecule  from  the 
species  N,  Ng  ,  or  NO  ,  can  be  used  to  represent  the  catalytic  body. 

If  M  is  the  only  such  fictitious  species,  then  N4  In  the  present 
example  Is  seven. 

The  conq)OSltlon  of  the  fluid  Is  specified  by  the  number  of  moles 
of  species  1  per  unit  mass  of  fluid  nj^  ,  and  is  referred  to  as  the 
mole-mass  ratio.  When  the  mole-mass  ratios  n^^  for  the  reacting 
species  (1  =  l,...,Ni)  are  known,  the  values  for  the  fictitious 
species  (i  =  Ni+  1,...,N4)  are  readily  found  from  alge¬ 

braic  equations.  In  addition,  there  exist  N^  -  N3  algebraic  equa¬ 
tions  that  express  conservation  of  components  and  the  constancy  of 
nj^  for  any  polyatomic  inert  species.  The  number  of  independent  values 
of  n^  ,  which  Is  the  minimum  number  that  must  be  accounted  for  by  rate 
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equationa.  Is  thus  -  (Nl-Nj^)  -  (Nj^-N3)  ■  Nj  .  The  additional  vari¬ 
ables  sure  given  by  algebraic  equations  of  the  form 


“l  “  ^lo"^  X  ^  “  N3+1,...,Nj^  ,  (U) 

>1 

where  the  aj^i  are  constants  and  the  Independent  species  are  numbered 
from  1  to  N3  .  The  first  Nj^-  Nj  of  equations  (4)  account  for  con¬ 
servation  of  components  and  Inert  species.  The  remaining  N^- 
equations  provide  for  any  nonslmple  catalytic  bodies. 

The  system  of  reactions  Is  represented  by 


Z  “rl*l 
1=1 


ri  1 


r 


where  represents  the  chemical  species  and  Qj.^  and  are  the 

stoichiometric  coefficients  of  the  reactants  and  products  respectively. 
The  forward  and  backward  rate  constants  for  reaction  r  are  denoted  by 
and  k^j.  • 

The  chemical  rate  equations  are  now  written  for  a  general  unsteady 
three-dimensional  flow  as  follows: 


C. 

I 

r=l 


ri  0  ’ 

r 


where  D(  )/Dt  is  the  substantial  derivative  and 


ri 


a 


ri 


1  -  1, 


A. 


(5) 


(6a) 


L^(p,T,n^)  ^  1  -  n  >  r  =  1,...,N2  ,  (6b) 


cr  k)«l 


[0r(p^T,n^)]‘^  s  (p)'*''  n  (6c) 


k=l 


8 


AEDC>TDR.«3.|2 


The  density  is  denoted  by  p  ,  the  equilibrium  constant  for  reaction  r 

by  K  _  ,  and  a  and  0  are  defined  as  follows: 
cr  r  r 


“r  '  I  “rj  '  Pr  ■  I  > 

J-1  >1 

For  convenience,  the  nuaber  of  rate  equations  (3)  1b  here  taken  as 
Instead  of  the  mlnlnum  nunber  N3  .  The  equations  expressing  conser¬ 
vation  of  components  will  be  taken  Into  account  later  when  the  number 
of  rate  equations  will  be  reduced  to  N3  .  We  may  also  note  that  emy 
catalytic  body  or  Inert  species  actually  appears  In  the  foregoing 
equations  only  In  0^.  ,  since  ■  0  for  1  «  1,...,Nij  .  Con¬ 

sequently,  all  Inert  species  are  taken  In  the  group  numbered  from 
Ni  +1  to  N]^  ,  formerly  used  only  for  fictitious  species.  The  first 
species  may  then  be  arranged  In  any  arbitrary  manner. 

3.2  LINEARIZED  FORM  OF  THE  CHEMICAL  RATE  EQUATIONS 

For  simplicity,  ve  require  that  the  reactions  be  sequenced  such 
that  the  first  are  linearly  Independent.  Using  the  notation 
(  )*  to  denote  tM  equilibrium  reference  state,  the  reference  compo¬ 
sition  is  then  given  by 

L^(p*,T*,n*)  =  0  ,  r  =  1,...,N3  ,  (8a) 


r  ■  1,. . . ,Ng  .  (7) 


and 


io 


I 

J-1 


ij  J 


i  =  N3+1,...,N^  ,  (8b) 


where  T  denotes  the  teiiq>erature ,  and  where  equations  (8b)  express 
conservation  of  components  (cf.  equations  (U)).  In  general,  as  noted 
eeurller,  the  quantities  with  asterisks  are  functions  of  position  and 
time.  These  quantities  are  thus  assumed  to  be  known.  Since  the  local 
equilibrium  conq>osltion  Is  also  necessary  for  the  linearization,  the 
pertinent  equations  cure  now  given  as  follows: 

I'j,(p>T,nj^^^)  =  0  ,  r  =  1,...,N3  ,  (9b) 

and 
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o 

"  ®lo  *  Z  ^  “  N3+1,...,N^  .  (9»)) 


>1 


Ve  now  perturb  about  the  equilibrium  reference  state,  denoting  the 
perturbation  quantities  by  (  )'  .  Thus,  ve  have 


p  =  p*  +  p*  , 

T  »  T»  +  T'  , 


n^  =  n*  +  n|  , 


i,e 


s  n* 


n*  +  n 


l,e  ' 


0  =  0*  +  6'  , 
r  r  r  ' 


1=  !,•••, , 


1=  , 


r  - 


(10) 


where  =  0  (p*,T*,n*)  .  Both  L  (p,T,nj^)  and  L  (p,T,n^  are  now 
expanded  about  the  reference  state thereby  resulting  In  * 

.  L^{p*,T»,n*)  P’ 

S/aL  \* 

*  Li4  . . 


and 


.  L^(p*,T*,n»)  P' 


I(^)  “i,. 


^  •  •  • 


r  -  l,...,ll2  .  (12) 


According  to  equation  (9a),  the  left-hand  side  of  equation  (12)  Is  zero. 
Subtracting  equation  (12)  from  (ll)  thereby  results  In  the  following: 
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(“i-  “i,.) 


T  ■  1^  •  a  •  f  aOlS®) 


Using  equations  (10) 


“i  ■  “i,e  -  -  “l,e  > 

ve  obtain  finally  the  llneaor  approoclmatlon 


H 


i  ■  l>a  a  a  f 


'  ■  l(%)  ' 


r  -  l,aa.,l^  a  (iSt) 


The  right-hand  side  of  equations  (5)  can  thus  be  linearised,  resulting 
In 


Dn, 

Dt 


1  -  1,...,N^  a  (15a) 


By  Interchanging  the  order  of  the  r  emd  i  susnatlons,  we  can  write 
equation  (l5a)  as  follows 


Ml 


1  ■  l,aaa,l^  f  (15b) 


where 


‘li 


«1|(P*,T*) 


■  I  ^(3^  ' 

x^l 


(16) 


Conservation  of  cooqponents  Is  finally  accounted  for  by  replacing  both 
n^  ^  and  n^^i  ■  Nj-i-1, . . .  ^R^)  In  the  right-hand  side  of  equation  (l5b) 
by'equatlons  (k)  ai^  (96).  After  some  rearrangement,  the  following 
Is  obtained: 
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1  ■  ,  (17) 


vhere 

h 

*11  ■  '*li 

>113+1 


1  ■  1^>..,H3^  i  ■  l^...fll3  .(18) 


3.3  TRANSFORMATION  OF  THE  CHEMICAL  RATE  EQUATIONS 

In  order  to  transform  the  linearized  form  of  the  chemical  rate 
equations  (I7)  to  their  canonlced  form,  ve  hypothesize  the  linear 
transformation 


1*3 

q^^  ■  ^  B^j^njj  ,  1  ■  1,...,N2  ,  (l9a) 

Ibil 

cmd  hence 

«3 

Ql^e  “  X  ®lk“k,e  *  1  -  l^.-.^Nj  ,  (19b) 

k»l 


vhere  the  q^  sure  generalized  mole-mass  ratios  and  the  sure 

constsuits.  ^e  canonical  form  of  equations  (17)  corresponding  to  equa¬ 
tions  (1),  is  then  given  by 

Dqi 

dT  “  ‘  '*i^‘^l"^i,e^  *  ^■'•••>1*3  >  (20) 

where  the  icj  are  functions  of  p*  smd  1*  .  So  fsur,  the  and 

sure  still  undetermined.  They  are  found  by  the  following  procedure: 

(a)  Differentiate  equation  (I9ci),  and  replace  Dnj^/Dt  by  mesms 
of  equation  (17),  thereby  obtaining 
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1>1  Ihl 


i  -  I,...,!,  .(21a) 


1%en  Intercbange  the  i  unBi  k  sumnatlona  to  obtain 


^?lkj  <“i,e“ 


n^) 


1  ■  l^...jNj  •  (21b) 


(b)  Replace  **  ^  in  equation  (20)  by  neans  of  equations  (19)  > 

thereby  olnalnlAg 


k-l 


l-l,...,Mj.  (22a) 


Replace  the  k-veurlable  of  sunnatlon  by  I  In  equation  (22a) 
as  follows: 


1-1 


1  -  1,...,N5  .  (22b) 


(c)  Compare  equations  (21b)  and  (22b)  to  obtain 


or 


1  -  1  -  1,...,N3  ,  (23) 
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where  le  the  Kroneeker  delta. 

(d)  For  the  monent,  aeauae  that  the  Kj^  and  K]^j  have  known 
eonetant  values.  With  the  i  eubseript  In  equations  (25) 
held  fixed,  we  then  have  Rj  homogeneous  linear  equatlcms 
for  the  R3  unknown  .  This  system  has  a  non-trlvlal 
solution  If,  and  only  If,  the  determinant  of  the  coefficients 
satisfies  the  condition 

*  0  >  1  •  1>*»*»J^3  •  (2^) 

These  N3  equations  cure  equivalent  to  the  single  character¬ 
istic  eqtiatlon 

-  0  ,  (24b) 

whose  N3  roots  are  the  Nj  eigenvalues  .  In  general, 
the  are  functions  of  p*  cmd  T*  and  therefore  can  be 

determined  once  an  equilibrium-flow  solution  Is  available. 

For  a  given  value  of  p*  cuid  T*  ,  the  nunibers  iCy ,  then 
determine  the  eigenvalues  by  means  of  equation  (24b}.  The 
,  In  conjunction  with  a  specific  ,  for  exeuiq>le  , 
are  next  used  to  determine  the  coefficients  (ki=l,. . .  ,N3) 
via  equations  (23).  Thus,  each  eigenvalue  Is  associated 
with  an  eigenvector  B^j^  .  The  resxxltlng  transformation  Is 
valid  however,  only  In  a  region  of  the  flow  that  Includes  the 
specific  values  chosen  for  p*  and  T*  .  The  extent  of  this 
region  depends  on  the  magnitude  of  the  changes  In  the  refer¬ 
ence  state.  For  example.  If  the  reference  state  Is  taken  to 
be  the  uniform  undisturbed  equilibrium  condition  for  a 
slightly  vcu:^lng  flow,  then  the  trcmsformatlon  Is  valid  for 
the  entire  flow  field. 

(e)  For  a  given  root  of  equation  (24b),  equations  (23)  do  not 

uniquely  determine  the  .  This  Indeterminacy  also  occurs 

In  vibration  theory.  One  possible  slnqple  way  of  removing  It 
Is  to  require  that 

>11  ■  1  i  1 . "5  • 


Thus  we  have  proved  the  existence  of  the  linear  transformation  (l^n)  that 
transforms  the  chemical  rate  equations  Into  their  canonical  form  (20). 

Two  points  still  remain  to  be  discussed.  First,  are  the  real 
and  positive?  This  Is  answered  In  the  affirmative  by  the  proof  given 
In  Appendix  I.  The  physical  significance  of  this  result  Is  discussed  In 
the  next  chapter. 
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Ibc  aeeond  point  is  that  ttas  trsnsforastlon  (I9<t)  can  be  extended 
to  the  aole-aass  ratios  n.  (i  *  N.-fl,. . . ,H.)  by  defining  new  gener¬ 
alised  mole -Bass  ratios  as  rollows: 


>  1  ■  Nj+  1, . . .  .  (26) 


Equations  expressing  conservation  of  con^nents  can  then  be  written  in 
terms  of  the  generalised  mole-mass  ratios  providing  no  two  eigenvalues 
are  equal,  as  is  generally  the  case.  (See  Goldstein  (1950)  for  addition¬ 
al  details.)  The  transfoimation  (l9s)  is  then  non-singular  and  an 
inverse  transformation  exists.  When  the  Inverses  of  transformation  (l9a) 
and  (26)  are  substituted  into  equations  (4),  the  desired  result  Is 
obtained  in  the  form 


qi  -  ^  A^^qj  ,  1  -  V  »  ^27) 

>1 

where  the  are  constants. 

4.0  GENERAL  PROPERTIES  OP  STIFF  EQUATIONS 

In  this  chapter,  certain  general  properties  of  stiff  equations  eure 
discussed.  These  properties  are  based  on  equations  (2)  and  cure  exten¬ 
sively  used  in  the  rest  of  the  report. 

For  convenience,  we  now  combine  equations  (2a)  and  (2b)  to  give 

(=8) 

where  e  is  a  small  constant.  The  psurameter  c  thus  corresponds  to 
UT]^  of  a  vibrational  rate  equation  or  to  u/k^  of  a  transformed  chem¬ 
ical  fate  equation.  The  general  solution  of  eouatlon  (28),  In  non- 
dimensional  form.  Is  readily  seen  to  be 

wtere  y(0)  ■  y^^^  Is  the  Initial  condition.  Sketch  1  shows  y/c  versus 
x/e  for  positive  c  for  two  different  values  of  yi^/c  .  The  line 
(y/®)  ■  ■  (x/®)  t  referred  to  as  the  local  equilibrium  curve,  corresponds 
to  qjj^  “  *11,0  *  Furthermore,  the  local  equilibrium  curve  corresponds  to 
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the  equlllbrluB-flow  solution  when  c  -*>0  .  The  other  dashed  curve « 

(y/c)  ■  -  (x/c)  +1  ,  referred  to  as  the  asymptote,  is  a  particular 
solution  of  equation  (28)  for  (yio/c)  ■  1  .  This  curve  corresponds 
to  the  true  nonequlllhrlua  solution.  Also  shown  In  Sketch  1,  hy  Iwavy 
line  segnents.  Is  the  direction  field  of  equation  (28)  for  three 
different  values  of  y/c  at  (x/c)  >  1  .  Because  of  the  converging 
nature  In  the  positive  x-dlreetlon  of  the  direction  field,  all  solu¬ 
tions  tend  toward  the  asymptote  exponentially  fast,  regardless  of  the 
Initial  value  yj^Q  .  As  a  result,  the  two  solutions  shown  differ  only 
sllcditly  fron  the  asymptote  for  any  x/c  greater  than  about  3  .  Any 
solution  thus  consists  of  two  portions:  a  transient  region  In  which 
the  solution  decays  rapidly  towards  the  asymptote,  and  a  subsequent 
region  in  which  the  solution  differs  only  negligibly  from  the  asymptote. 
When  equation  (29)  Is  plotted  in  dimensional  coordinates  (i.e.,  x  and 
y  )  and  c  Is  very  small  (e.g. ,  e  ■>  10 ),  then  the  local  equilibrium 
curve  y  *  -  X  lies  very  close  to  the  asyiiq>tote.  Consequently,  the 
direction  field,  when  e  Is  small,  changes  extremely  rapidly  In  th* 
vicinity  of  the  local  equilibrium  curve. 

Sketch  2  Is  similar  to  Sketch  1  except  e  Is  now  negative.  The 
coordinates  In  this  sketch  are  chosen  as  -  (y/c)  cuid  -  (x/c)  In 
order  that  the  curves  (y/c)  ■  -  (x/c)  aM  (y/c)  =  -  (x/c)  +  1  have 
similar  orientation  In  the  two  figures.  The  direction  field  now  diverges 
In  the  positive  x-dlrection,  as  shown  at  -  (x/c)  *  1  .  Thus,  eLLl 
solutions  tend  to  diverge  exponentially  fast  from  (y/c)  ■  -  (x/c)  +  1  , 
which,  for  convenience,  is  still  referred  to  as  the  asymptote. 

From  a  physical  point  of  view,  negative  values  for  c  (or  ) 
Inqply  that  the  nonequillbrlum-flow  solution  cannot  remain  close  to  the 
equlllbrluffl-flow  solution.  Any  nonequilibrium  nozale  cailculation,  even 
with  equilibrium  Initial  conditions,  would  In  this  case,  diverge  Ismiedl- 
ately  from  the  equilibrium-flow  solution.  This  situation  Is,  however, 
physically  unrealistic.  We  have,  therefore,  established  the  fact  that 
negative  c  (or  k.  )  cannot  occur  by  physical  reasoning.  Appendix  I 
also  establishes  this  by  purely  mathematlced  reasoning.  Consequently, 
positive  c  will  be  assumed  In  the  rest  of  this  report. 


5.0  RUNCE'KUTTA  ANALYSIS 

The  Runge-Kutta  procedure  Is  one  of  the  oldest  and  most  commonly 
used  methods  for  the  numerical  integration  of  ordinary  differential 
equations.  This  chapter  considers  the  behavior  of  this  procedure  as 
applied  to  the  Integration  of  the  stiff  equation  (28).  Our  main  result 
will  be  the  determination  of  the  truncation  error  and  of  a  stability 
criterion.  To  clarify  the  discussion,  a  brief  preliminary  description 
of  the  Runge-Kutta  procedure  will  be  helpful.  For  additional  details, 
the  Interested  reader  is  referred  to  Ince  (1926),  Hildebrand  (1956)  or 
Ralston  emd  Wilf  (i960). 
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The  Runge-Kutta  procedure  Is  designed  to  approximate  the  Taylor- 
serles  solution 


y(Xp+h)  -  y(x^)  +  hy'(x^)  +  —  y"(x^)  +  ...  , 


(30) 


to  the  equation 


y'(*)  =  ^  “  f(x»y)  > 


(31) 


at  the  point  Xq+  h  .  The  qucmtlty  h  denotes  the  integration  step  size. 
In  contrast  to  the  formal  Taylor-serles  solution  (30),  the  Runge-Kutta 
procedure  does  not  require  the  explicit  evaluation  of  the  second  or 
higher  derivatives.  Rather,  an  approximation  to  expansion  (50)  Is  obtain¬ 
ed  at  the  expense  of  several  evaluations  of  the  first  derivative.  The 
approximation  to  the  first  p -t-l  terns  on  the  rl£^t-hand  side  of  (30), 
referred  to  as  a  p^^-order  Rvuige-Kutta  procedure,  requires  the  evaluation 
of  p  first  derivatives.  These  derivatives  are  evaluated  at  various 
points  In  the  x,y  plane  as  follows: 

^2  ’  ‘■'V  V,  V  ' 

*■3  •  '<V  “3“'  V  P3l“'l  +  ®32“2>  ’  )  (32) 


f  =  f (x  +  ah,  y  +h  >  p  ,f .)  , 
p  o  p  "  ■'o  L.  Vi  y  ’ 

J-1 


where  y^  ^  y(x-)  and  the  and  are  constants  whose  valties  are 

discussed  shortly.  The  desired  approxlmtlon  is  then  obtained  from 


y(Xo+  h) 


=*  ^  L 


(33) 


1=1 


The  constants  7^  >  end  p^j  are  determined  by  the  condition  that 
equations  (30)  and  (33)  match  to  the  p^^  order  at  the  point  (xp,yQ)  . 

In  other  words,  the  coefficients  of  like  powers  of  h  from  (hjo  to 
(h)P  are  equated  with  each  other.  This  matching  condition  does  not 
determine  all  the  constants.  For  exEuqple,  for  a  fourth -order  procedure. 
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the  matching  condition  leaves  two  of  the  constants  vinspecified.  Values 
for  these  two  constants  nay  therefore  be  choseh  arbitrarily. 

At  the  start  of  a  nonequilibrium  calculation,  as  Section  ^.1  will 
demonstrate,  the  step  else  h  must  be  smedl.  It  is  therefore  Inqperative 
that  provision  be  Included  in  the  integration  procedure  for  increasing 
or,  if  necessary,  decreasing  the  step  size  as  the  calculation  proceeds 
in  the  X'dlrection.  Any  technique  for  vaucylng  h  requires  an  estimate 
of  the  error  that  is  introduced  into  the  solution  when  a  single  step  is 
integrated.  This  error  is  actually  an  estimate  of  the  combined  round¬ 
off  and  truncation  errors.  To  estimate  this  error  for  a  step  of  size 
h  ,  the  step  is  reconputed  by  means  of  tvo  steps,  each  of  length  h/2  , 
and  the  results  conpared. (See  Ince  (19S6)  for  further  details.)  Althou^ 
this  is  e\istomarily  referred  to  as  a  single  step  in  the  integration,  the 
numerical  solution.  In  fact,  is  primarily  given  by  the  tvo  half^step  cal¬ 
culations,  vhlch  are  more  accurate  than  the  full-step  result.  (This 
terminology  is  also  used  in  Part  II.)  Each  step  of  length  h  thus 
requires  3p  complete  evaluations  of  the  right-hand  side  of  the  differ¬ 
ential  equations  if  the  error  estimation  is  performed.  In  compeurison 
with  tvo  half  steps  vlthout  error  estimation,  vhlch  vould  require  2p 
evaluations,  ve  see  that  error  estimation  increases  overall  computation 
time  approximately 

Runge-Kutta  procedures  vlll  be  conpau*ed  vlth  predictor-corrector 
procedxires  in  Chapter  7*0.  These  latter  procedures,  hovever,  do  not 
estimate  the  error  by  recomputing  vlth  tvo  half  steps.  This  comparison 
therefore  requires  that  each  half  step,  henceforth  denoted  by  h  ,  be 
considered  as  the  Runge-Kutta  step  size*.  The  number  of  eveduatlons  of 
the  right-hand  side  of  the  differential  equations,  in  this  Instamce,  is 
taken  as  l.^P  per  integration  step. 

In  practice,  the  most  vldely  used  Runge-Kutta  procedure  is  the 
fourth -order  one.  The  results  obtained  from  the  subsequent  analysis 
vlll  be  valid,  hovever,  for  any  order  Runge-Kutta  procedure  from 
second-order  (  p  «  2  )  on.  Because  of  the  considerable  complexity 
of  this  procedure,  such  generality  is  rare.  For  an  exanple  of  hov 
extremely  complicated  the  results  usually  are,"  the  reader  is  referred 
to  a  paper  by  Lotkin  (1951).  The  error  associated  vith  the  inte¬ 
gration  of  equation  (31)  by  a  specific  fourth-order  Runge-Kutta  pro¬ 
cedure  is  given  in  that  paper.  The  simplicity  and  generality  of  the 
results  that  will  be  obtained  here  are  directly  due  to  the  simple 
mathematical  structure  of  equation  (28). 


* 

A  single  step  of  length  h'  in  Part  II  vlll  now  be  considered  as  tvo 
steps,  each  of  length  h  ,  vhere  h'  >  2h  . 
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S.1  TRUNCATION  IRROR 

We  nov  ehov  that  the  values  chosen  for  the  free  constants  in 
the  Runge-Kutta  procedure  are  Innaterlal  when  f  Is  given  by* 


f  - 


.  y-*-* 

c 


(3U) 


The  first  step  in  the  proof  Is  to  show  that  equation  (33)  does  not 
contain  higher  powers  of  h  than  (h)^  .  To  do  this,  substitute  equa¬ 
tion  (3h)  Into  (32)  as  follows: 


f^  -  -  €-^(y^+  x^) 
f 2  =  “  *  (yo'*'  ^21^1^  *0'*’ 


V  V 


(55) 


Thus,  does  not  Involve  h  ,  f2  Is  linear  In  h  ,  f3  is  quadratic 
In  h  ,  etc.  The  highest  power  of  h  in  fp  is  therefore  (h)P“^  . 

Since  the  summation  term  In  equation  (33)  Is  multiplied  by  em  h  ,  the 
highest  power  of  h  In  this  equation  is  (h)P  .  Consequently,  the  rlght- 
hemd  side  of  equation  (33)  is  precisely  equal  to  the  sum  of  the  first 
p+1  terms  of  the  Taylor-serles  expansion  (30).  Thus,  no  matter  what 
values  are  chosen  for  the  free  constants,  equation  (33)  will  compute  the 
same  number,  thereby  proving  the  above  assertion.  The  underlying 
reason  for  this  result  is  that  f  ,  as  given  by  equation  (34),  is  linear 
In  both  X  and  y  . 


As  a  consequence  of  the  foregoing,  a  simple  expression  for  the 
truncation  error**  can  now  be  obtained.  The  truncation  error  is  defined 
as  the  difference  between  the  (p+2)^^  term  in  the  Taylor-serles 
expansion  (30)  and  the  corresponding  term  coopted  by  equation  (33). 

The  (p  +  2)^^  term  in  the  series  expansion  is  given  by 


*  /  X  -1 

The  quantity  -(e)  ,  as  used  in  this  report,  is  equivalent  to 

A  5  (df/dy)  ,  which  Is  frequently  used  in  the  literature  of  numerical 

analysis. 

**  This  truncation  error  should  not  be  confused  with  the  estimate  for 
the  combined  round-off  and  tzouicatlon  errors  mentioned  In  Section  ^.0. 
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,4xr*V 


whereas  the  corresponding  term  of  equation  (33),  l.e.,  the  term  contain¬ 
ing  (h)^'*'^  ,  is  zero.  Denoting  the  truncation  error  of  a  p'*'^-order 
Runge-Kutta  procedure  by  Ep  ,  ve  therefore  obtain  the  result 


7(^  >V  -o  -  (-  ' 


where  equation  (34)  is  used  to  evaluate  bhe  derivative  in  expression  (36). 

Three  Important  conclusions  stem  from  this  equation.  First,  the 
term  Yq*  Xq  ~  e  represents  the  vertical  distance  from  a  point 
(xq,-  Xq  +  c)  on  the  asymptote  to  the  integration  point  (xo,yQ)  . 

Thus,  when  the  numerical  solution  is  close  to  the  asymptote  this  term 
is  small.  In  this  case,  it  is  possible  to  have  h  >  e  and  still  have 
a  small  truncation  error.  A  necessary  condition  then  for  a  large  inte¬ 
gration  step  size  is  that  the  numerical  solution  be  close  to  the  asymptote. 
In  terms  of  the  integration  of  the  nonequilibrium  equations  of  Part  II, 
this  means  that  the  numerical  solution  must  be  close  to  the  physical  solu¬ 
tion  in  order  that  the  step  size  be  large. 

Second,  when  Xq  »  y  ■  0  the  tr\incation  error  is  given  by 


The  condition  Xq  >  y^  ■  0  represents  the  equilibrium  initial  condition 
for  the  rate  equations.  This  topic  was  discussed  in  detail  in  Part  II. 

We  here  note  that  for  small  Ep  ,  in  this  condition,  it  is  necessary  that 
h  be  less  than  e  .  For  example,  with  p  >  4  and  |Ep|  «  10~^  ,  equa¬ 
tion  (38)  requires  that  (h/s)  ■  .2605  • 

The  third  result  to  stem  from  equation  (37)  is  that  the  condition 


requires  that 
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h/€  S  p  +  2  .  (40) 

In  other  words,  for  the  same  integration  step  size,  higher-order  proce¬ 
dures  will  have  smaller  absolute  truncation  errors  provided  condition 
(UO)  is  satisfied.  This  condition  must  be  fulfilled,  as  will  be  shown, 

In  order  that  the  Integration  procedure  remain  stable. 

S.2  STABILITY 

This  chapter  concludes  with  a  discussion  of  stability  of  Runge- 
Kutta  procedures  applicable  to  the  region  in  which  the  n\unerlcal  solution 
differs  only  negligibly  from  the  asymptote. 

The  accuracy  of  a  numerical  integration  depends  both  on  the  error 
incurred  at  each  step  of  the  Integration  and  on  how  the  error  grows. 

If  the  error  decreases  as  the  succeeding  steps  are  performed,  the  Inte¬ 
gration  procedure  is  referred  to  as  stable.  In  other  words,  If  the 
numerical  solution  converges  to  the  exact  solution  as  the  succeeding 
steps  are  performed,  the  Integration  procedure  Is  considered  to  be  stable. 
If  the  error  remains  constant,  the  procedure  is  termed  neutrally  stable, 
and  If  the  error  increases,  the  procedure  is  termed  unstable. 

We  now  turn  to  the  problem  of  applying  this  definition  of  stability 
to  our  consideration  of  equation  (28).  As  Sketch  1  shows,  all  analyt¬ 
ical  solutions  tend  toward  the  asymptote.  Consequently,  any  stable  nu¬ 
merical  solution  must  also  tend  toward  It.  In  other  words,  the  distance 
from  the  point  (Xq+  h,y(xQ+  h))  to  the  asymptote  must  be  less  than 
the  distance  from  (xQ,yQ)  to  the  asymptote,  in  order  for  the  numerical 
procedure  to  be  stable.  For  simplicity  of  formulation,  we  shall  use 
vertical  distance  rather  than  actual  distance;  the  final  result,  of 
course,  is  the  same.  Sketch  3  presents  a  sketch  of  the  pertinent 
quantities.  The  two  heavy  vertical  lines  represent  the  two  distsmces. 
Their  ratio,  denoted  by  6  ,  is  given  by 


-(x  +  h)+  €  -  y(x  +  h) 
o  'o 

where  the  numerator  Is  the  vertical  distance  at  Xo+  h  .  After  a  certain 
amount  of  algebraic  manipulation  6  can  be  written  as 


6 


1  - 


e  e 


1  - 


y  +  X 
'o  o 


(4lb) 
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The  quantity  y(xo+  h)-  y^  In  equation  (4lb)  Is  to  be  computed  by  a 
p^^-order  Runge-Kutta  procedure.  As  shown  previously,  the  same  value 
can  also  be  obtained  from  the  Taylor-serles  exi>anslon  (30)  truncated 
after  the  (h)P  term.  When  this  expansion  Is  written  out,  we  obtain 
the  equation 
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where  the  curly  brackets  contain  the  derivatives  /'(xq), 

*  evaluated  in  accordance  with  equation  (28). 
Again  after  some  algebraic  manipulation,  equation  (42a)  can  be  re¬ 
written  in  the  more  convenient  form 


'  '  1^1 


(42b) 


After  substituting  eqviatlon  (42b)  into  equation  (4lb),  we  obtain  the 
final  result 


-  1 


e  e  e 


(4lc) 


According  to  the  foregoing  considerations,  a  Runge-Kutta  procedure 
of  order  p  Is  stable  In  the  integration  of  equation  (28)  If 

l6|  <  1  .  (43) 

For  each  Integer  value  of  p  there  will  exist  two  non-negative  values 
of  h/c  such  that  |b|  ■  1  .  One  of  these  values  Is  (h/e)  =  0  . 

The  other  value,  denoted  by  the  stability  parameter  B(p)  ,  represents 
the  upper  limit  for  values  of  h/c  for  which  the  Integration  Is  stable. 
Thus,  condition  (43)  Is  equivalent  to  the  following  stability  condition 

h/c  <  B(p)  ,  (44) 

where  the  quantity  B(p)  is  determined  by  the  relation  |&(B;p)|  >  1  . 
The  following  table  lists  B(p)  for  different  values  of  p  . 
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STAilLITY  PARAMITER  B(p) 


p 

&(B;p) 

B(p) 

2 

1 

2.00 

3 

-1 

2.52 

k 

1 

2.00 

5 

-1 

3.20 

An  Important  conclusion  of  the  above  analysis  Is  that  stability 
considerations  Impose  a  limit  on  the  Runge-Kutta  step  size  when  c 
Is  small.  Stable  Integration  with  (h/e)  >  1  >  however.  Is  possible. 
When  (h/c)  >  1  ,  the  truncation  error  will  be  Inordinately  large 
unless  yQ+  Xq  -  e  Is  small,  l.e.,  unless  the  numerical  solution  Is 
quite  close  to  the  asymptote.  Thus,  only  after  the  transient  region 
Is  completed  are  large  Integration  steps,  l.e.,  B>h/€>1,  possible. 

Negative  values  for  B(B;p)  In  the  above  table  signify  that 
y(xo+  h)  and  Yq  are  on  opposite  sides  of  the  asymptote.  For  example, 
for  p  =  4  the  integration  procedure  is  stable  for  0  <  h/e  <  2.80  , 
and  yQ  and  y(xQ+  h)  are  both  on  the  same  side  of  the  asymptote. 

The  above  table  shows  that  condition  (40),  as  already  mentioned,  is 
always  satisfied  when  the  integration  procedure  is  stable.  Therefore, 

If  h/e  is  chosen  to  achieve  stability,  then  the  hi^er-order  the 
procedure,  the  smaller  the  truncation  error.  This  table  also  demon¬ 
strates  that  higher-order  Runge-Kutta  procedures  are  more  stable  them 
lower-order  procedures.  Each  Increase  In  order,  however,  requires 
additional  evaluations  of  the  derivative  f(x,y)  (see  equation  (31)). 
Since  the  evaluations  are  frequently  very  time  consuming  on  a  digital 
computer,  a  compromise  between  stability  and  computation  time  becomes 
necessary.  When  the  truncation  error,  stability,  cmd  computation  time 
cure  all  considered,  the  fourth-order  procedure  appears  to  be  a  reason¬ 
able  choice. 

The  value  of  B  ,  given  in  Part  II  for  a  fourth-order  Runge-Kutta 
procedure.  Is  based  on  a  double  step  of  length  2h  .  Its  value  Is 
therefore  2  x  2.0  ■  5.6  .  (The  value  actually  given  in  Part  II  is 
5.7;  and  Is  slightly  in  error.) 
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4.0  KIDICTM^MRECTOR  ANALYSIS 

This  chapter  analyzes  predictor -corrector  procedures  In  a  feushlon 
dlmllar  to  that  of  the  Runge-Kutta  emalysls.  As  In  the  preceding 
chapter,  a  brief  description  of  predictor -corrector  procedures  will 
help  clarify  the  analysis.  For  additional  details,  the  reader  Is 
referred  to  Ralston  and  Vllf  (i960)  or  Heusnlng  (1962).  The  latter 
reference,  In  psurtlcular,  contains  an  extensive  discussion  of 
predictor-corrector  procedures. 

Any  predictor -corrector  procedure  Is  a  finite-difference  method 
that  uses  previously  conqinited  data  to  assist  In  obtaining  yCx^-t-  h)  . 

In  other  words,  the  calculation  of  yCx^-t-  h)  generally  requires  the 
use  of  y(xQ),  y(xQ-  h),  y(xQ-  2h),  ...  and  of  the  derivatives 
y'(xQ+  h),  y'(xQ),  y'(xQ-  h),  ...  .  The  more  back  data  used,  the 

higher  Is  the  order  of  the  procedure.  This  use  of  back  data  repre¬ 
sents  the  chief  difference  between  a  Runge-Kutta  procedure,  which  Is 
not  a  finite-difference  method,  smd  a  predictor-corrector  procedure. 
Nevertheless,  both  types  of  procedure  attempt  to  estloiate  the  solution 
of  a  differential  equation  by  means  of  a  Taylor-serles  expansion. 

A  typical  corrector  (or  predictor)  formula  Is  given  by  (see 
Hamming  (1962)) 

y(x^+  h)  = 

+  hrb_^y'(x^+  h)  +  b^y'(x^)  +  b^y'(x^-  h)  +  b2y'(x^-  2h)|  , 

(45) 

where  a^  and  bj^  are  constants.  Of  the  seven  constants  here,  five 
are  usually  chosen  such  that  formula  (U5)  matches  the  Taylor-serles 
expansion  (30)  to  fourth-order.  The  remaining  two  constants  may  then 
be  chosen  to  improve  stability.  Corrector  formulas  require  the  use  of 
y'(xQ+  h)  ,  l.e.,  b_]^  4  ^  >  whereas  predictor  formulas  do  not,  l.e., 

b_l  =  0  .  With  b_j^  =  0  /  equation  (45)  Is  first  used  to  predict  a 
value  for  y(xQ-i-  h)  .  Based  on  this  predicted  value,  a  first  estimate 
for  the  derivative  y'(xQ+  h)  is  computed  by  means  of  the  differential 
equation.  This  derivative  Is  next  used  to  calculate  a  new,  more  precise 
value  for  y(xQ+  h)  by  means  of  a  corrector  formula.  A  second,  EUid 
final,  value  for  y'(xQ+  h)  is  now  computed  using  this  new  value  for 
y(xQ+  h)  In  conjunction  with  the  differential  equation.  The  last  values 
for  y(xo+  h)  and  y'(xQ+  h)  are  the  ones  used  for  the  calculation  of 
the  next  step. 

As  already  noted  In  Chapter  5*0,  any  feasible  Integration  procedure 
must  be  able  to  vary  the  integration  step  size.  Therefore,  a  method 
for  estimating  the  local  truncation  error  Is  necessEoy.  An  error  estimate 
with  the  above  predictor-corrector  procedure  Is  possible  If  the  order 
of  the  truncation  errors  for  the  predictor  and  the  corrector  are  equal. 
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We  shall  always  assxjne  this  to  be  the  case.  The  truncation  error  of 
the  overall  procedure  is  then  of  the  same  order  as  that  of  the  predictor 
or  corrector.  This  error  Is  obtained  by  properly  con^rlng  the  results 
of  the  predictor  and  the  corrector  calculations.  One  often  quoted 
advantage  of  predictor-corrector  as  conqMred  to  Runge-Kutta  procedures 
is  that  no  additional  evaluations  of  the  derivative  f(x,y)  (see 
equation  (31))  are  necessary  to  estimate  the  error. 

The  analysis  contained  In  this  chapter  Is  restricted  to  the  fore¬ 
going  predictor-corrector  procedure.  Unlike  Runge-Kutta  procedures, 
there  are  many  different  types  of  predictor-corrector  procedures.  An 
exhaustive  analysis  of  the  numerous  types  Is  beyond  the  scope  of  this 
work.  In  view  of  this,  the  results  of  this  chapter  are  of  a  preliminary 
nature  only. 

A  further  comparison  between  Runge-Kutta  and  predictor-corrector 
procedures  will  be  useful  before  discussing  truncation  error  and  stability. 
First,  predictor-corrector  procedures  are  not  self -starting.  Thus,  a 
Runge-Kutta  procedure,  which  Is  self-steurtlng.  Is  generally  used  to  com¬ 
pute  the  first  few  steps.  As  a  consequence,  predictor-corrector  pro¬ 
cedures  are  more  Involved  than  Runge-Kutta  procedures  to  program  for  a 
computer. 

Second,  varying  the  Integration  step  size  is  more  troublesome  with 
predictor -collector  procedures,  where  the  step  size  Is  either  halved  or 
doubled,  as  compared  to  Runge-Kutta  procedures,  where  the  step  size  can 
be  changed  arbitrarily.  According  to  criterion  (4h),  which  is  also 
applicable  here,  a  predictor-corrector  procedure  will  be  stable  after 
the  step  size  Is  doubled  only  If  (h/s)  <  ,  where  h  Is  the  orig¬ 

inal  step  size.  Thus,  the  step  size  must  be  less  than  half  Its  permis¬ 
sible  size  If  the  procedure  is  to  be  stable  after  the  step  size  is 
doubled.  The  stability  parameters  that  will  be  derived  in  Section  6.2 
therefore  overestimate  the  average  step  size  that  Is  possible.  An 
efficient  technique  for  controlling  the  step  size  would  double  it  when¬ 
ever  (h/e)  <  .  In  this  circumstance,  the  step  size  would  on  aver¬ 

age  be  (3/4)  that  allowed  by  criterion  (44).  Consequently,  a  more 
realistic  stability  criterion  for  predictor-corrector  procedures  Is 
given  by 


as  compeared  to  criterion  (44). 

Third,  the  nature  of  instability  Is  different  for  the  two  procedures. 
As  we  have  seen,  Runge-Kutta  Instability  is  due  to  the  nxunerical  solution 
diverging  from  the  exact  one.  Any  predictor-corrector  procedure,  being 
a  finite -difference  method,  may  converge  toward  one  of  perhaps  several 
spurious  solutions  instead  of  toweurd  the  exact  solution.  When  this 
occurs  the  procedure  is  also  said  to  be  unstable.  Both  types  of 
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instability,  however,  stem  from  too  large  a  value  for  the  parameter 
h/e  .  Thus,  stability  considerations  impose  a  limit  on  the  step  size 
for  both  types  of  integration  procedures. 

Finally,  predictor -corrector  procedures  require  fewer  evaluations 
of  the  derivative  f(x,y)  and,  therefore,  may  require  less  computation 
time  them  Runge-Kutta  procedures.  In  particular,  the  foregoing  pre¬ 
dictor-corrector  procedure  requires  only  two  evaluations  of  the  deriv¬ 
ative  per  integration  step.  With  emy  Runge-Kutta  procedure,  the  corre¬ 
sponding  number  of  evaluations  is  1.5p  • 

6.1  TRUNCATION  ERROR 

The  truncation  error  for  a  predictor-corrector  procedure  of  order 
p  Is  given  by 


E 


P 


c 


P 


(e) 


(47) 


where  Cp  Is  a  constant  that  depends  on  the  specific  procedure.  The 
derivative  y(P"*'^^(0)  Is  evaluated  at  a  point  6  in  the  interval 
from  Xq-  rh  ,  the  farthest  back-data  point  used  In  the  integration 
formula,  to  Xg+  h  .*  This  derivative  is  essentially  the  same  as  that 
in  expression  ([36)  with  the  exception  that  here  x  -  6  .  Thus,  when 
equation  (28)  is  Integrated,  the  truncation  error  is  still  given  by 
equation  (37),  except  for  the  constant  c-  ,  which  generally  is  within 
an  order  of  magnitude  of  unity,  and  the  replacement  of  yQ+  x^  by 
y(0)  +  9  ,  Consequently,  the  first  two  conclusions  in  the  previous 
section  following  equation  (37)  edso  appl>  here.  In  particular,  it  is 
again  possible  for  h  to  be  greater  than  e  and  for  the  truncation 
error  nevertheless  to  be  small  provided  the  numerical  solution  is 
sufficiently  close  to  the  asymptote. 

6.2  STABILITY 

At  least  two  different  definitions  of  stability  for  predictor- 
corrector  procedures  are  currently  in  use.  The  first  of  these,  referred 
to  as  relative  stability,  requires  that  spurious  solutions  be  small 
relative  to  the  exact  solution.  Relative  stability  is  important  only 
when  h/c  is  negative  and  therefore  is  not  considered  here.  This 
section  deals  with  the  more  common  form  of  stability,  which  is  not  as 
restrictive  as  relative  stability,  but  is  applicable  only  when  h/e 
is  positive.  A  finite -difference  scheme  is  considered  stable  when  all 
the  roots  of  its  characteristic  polynomial  are  interior  to  the  unit 
circle.  In  this  circumstance,  all  spurious  solutions  of  the  flnlte- 


*  Equation  (47)  assumes  that  the  influence  function  is  of  constant  sign 
in  the  interval  from  Xq-  rh  to  Xq+  h  .  For  additional  details  see 
Hamming  (1962). 
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difference  schene  decay  toward  zero.  Thus,  once  the  transient  region 
(see  Chapter  U.O)  Is  completed,  the  nunerlcal  solution  must  converge 
to  the  exact  solution  when  the  finite -difference  scheme  Is  stable. 


There  are  many  stability  studies  of  the  type  of  predictor-corrector 
procedure  that  we  are  here  concerned  with.  For  Instemce,  the  reader 
Is  referred  to  Dahl^ulst  (1936)  and  (1939)^  Hamming  (1939)  end  (1962), 
Crane  and  Lambert  (1962),  Henrlcl  (1^2),  and  Chase  (1962).  With  the 
exception  of  Chase,  all  of  the  referenced  authors  have  assumed  that 
stability  of  the  overall  predictor-corrector  procedure  Is  determined 
primarily  by  the  stability  of  the  corrector  alone.  Consequently,  these 
authors  derive  and  examine  the  characteristic  polynomial  relevant  only 
to  the  corrector.  The  foregoing  assumption,  however.  Is  usually  correct 
only  In  the  limit  as  h  tends  to  zero.  By  Investigating  a  few  special 
predictor-corrector  procedures,  Chase  (1962)  has  shown  that  this  assump¬ 
tion  may  not  be  valid  as  h/e  Increases.  Therefore,  to  determine  sta¬ 
bility  parameters,  it  Is  necesseury  to  consider  the  overall  integration 
procedure,  not  Just  the  corrector.  This  is  accomplished  by  deriving 
a  characteristic  polynomial  for  the  overall  procedure. 


We  now  define  a  generalized  predictor-corrector  procedure  by  the 
following  formulas: 

k'  k’ 

z(x^+  h)  «  ^  ^  Yj  ®i^'^*o" 

i»0  1*0 


h) 


i»0 


ih) +  hb_^z 


(V 


h)  +  h 


A. 

I 

i-0 


'^V 


ih)  ,  (49) 


where  Aj,  a.  and  bj  are  constants 
Integers.  Equation  (48)  is  the  predictor 
first  estimate  z(xq+  h)  for  y(xQ+  h)  . 
the  corrector .  formula  (49) ,  is  given  by 
eqxiation  (3l)). 


and  k'  and  k  are  positive 
and  is  used  to  determine  a 
The  quemtity  z'(xq+  h)  in 
f(xjj+  h,  z(xq+  h))  (see 


As  mentioned  earlier,  we  shall  require  that  each  of  the  above 
formulas  have  a  truncation  error  of  order  (h)?"*"^  .  This  condition  will 
then  determine  some  of  the  constants  Ai,  Bj^,  a^^  ,  and  bj^  .  The  formu¬ 
lation  of  this  condition  also  requires  that  the  estimates  z(xq+  h)  and 
z'(xo+  h)  be  replaced  by  y(xQ+  h)  and  y'(xQ+  h)  respectively. 

When  this  is  done,  formulas  (4o)  and  (49)  are  of  the  same  form  except 
that  B_2  ■  0  in  formula  (48).  To  satisfy  the  condition  that  the 
corrector  be  of  order  p  ,  we  require  that  formula  (49)  be  exact  for 
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the  functions  y  ■  .  In  other  words,  the  p  +  1  functions 

y  -  1  ,  ^ 


y  ■  X 


J 


(50) 


sure  substituted  successively  into  formula  (^9),  thereby  resulting  in  the 
system  of  equations 


K  K 

(x^+  h)“  -  y  a^(x^-  ih)"+  mh  y  b^( 
1-0  1^-1 


..  vm-1 
x^-  ih) 


* 


(51) 


Since  the  aj^'s  and  b^'s  ought  to  be  Independent  of  the  arbitrary 
parameters  Xq  and  h  ,  it  should  be  possible  to  replace  equations  (91) 
by  a  simpler,  alternate  system  not  containing  these  two  parameters. 

This  system,  given  without  proof  by  Dahlqulst  (1956),  is  derived  in  the 
following  paragraph.  Eqviation  (51)^  and  hence  equation  (59) >  can  also 
be  used  to  determine  the  constemts  in  the  predictor  providing  aj^,  bj^  , 
and  k  are  replaced  by  A^,  ,  emd  k'  respectively  and  »  0  . 

We  proceed  as  follows;  Apply  the  binomial  theorem  successively 
to  (xo+  h)®  ,  (xq-  Ih)®  ,  and  (xq-  Ih)®"^  ,  thereby  resulting  in  the 
following  equations: 


'v*-'"- '  '52) 

J*0 

where 

^  j  j  s  binomial  coefficient  »  ®(®*^)  •  •  •.(®~J'*'^)  ^  (53) 
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(V  > 

J-O 


Y  (-l)J(l)J(;)x;-V  , 

>0 


and 


m-1 


J'-O 

Replace  J '  Ijy  J  - 1  ,  and  use  the  identity 

t:i)  -= 


(55) 


(56) 


to  transform  equation  (55)  into 


(x^-  ih)“"^  ”  ■  i  X  •  ^57) 


Substitute  equations  (^2), 

J-i 

(5^)  and  (57)  into  equation  (5I)  and  then 

Interchange  the  J 

and  1 

summations,  thereby  resulting  In 

m 

m 

k 

’  I 

(.i)j(j)xrv  £ 

J-O 

J-O 

1-0 

m 

k 

-E 

J=0 

i--l 

n^O^X^vas^P  •  (5®) 
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Vhen  the  coefficients  of  like  powers  of  h  ,  or  Xq  ,  are  equated,  the 
desired  result  Is  obtained: 


Z  ^  Z  J  •  (59) 

1^0  1^-1 


Equations  ($9)  have  the  following  sltqple  form  vhen  explicitly  written 
out 

W  V***'"  \ 

a^+  2a^+...+  kaj^  -  (b_^  +  +. . .+  bj^  )  «  -  1  , 

ai+  4a2+.  ..+k^aj^  -  2(-h_^  +  b^+  2h^  +...+  kbj^  )  -  1  ,  \(60) 


+kPaj^.p((-l)P"\j^+bj^+2^"^>2+...+(k)^"^>j^^  -  (-1)**  . 


To  fxirther  Illustrate  the  Importance  of  Chase's  assertion,  we  shall 
analyze  the  Adams  method  In  detail.  This  will  be  accomplished  by  exam* 
ining  the  stability  of  the  predictor  alone,  of  the  corrector  alone,  and 
of  the  overed.1  procedure.  Furthermore,  in  order  to  parallel  the  hunge* 
Kutta  analysis  given  In  Chapter  ^.0,  this  analysis  will  also  consider 
all  orders  from  p  ■  2  to  p  ■  5  • 

We  now  give  without  proof  (see  Chase  (1962))  the  characteristic 
polynomials  needed  for  the  Investigation.  All  of  the  polynomials  are 
based  directly  on  formulas  (48)  and  (49).  The  various  constants  that 
appear  In  them  will  then  be  specialized  for  the  Adams  method  wd  for 
equations  (60).  The  polynomials  are  denoted  by  C'P*"),  and 

C^P®)  where  the  superscripts  denote  predictor,  corrector,  and  predictor- 
corrector  respectively.  They  are  as  follows: 


c(P^)  E  +  ^'(.A^+  |)  .  0  ,  (61) 

1=0 


32 


AeDC*TDR.«3.l3 


,(co) 


(1  +  ^  (-a^+  -  0  , 


(62) 


1^0 


1-0 


c'"”'  E  < 


0  1^0 

if  k  i  k’,  (63) 

k  k' 

V/  hvk’-l^.  ^VfA  TJ  hv  k'-i  - 

p  +2.^‘V^  ^ ^-1  I^P  *  °  * 

if  k  S  k'  .  (6U) 


1^0 


i«0 


Following  Nordsleck  (1962),  we  shall  require  any  Adanus  formula  to 
have  a  characteristic  i)Olynomlal  of  the  form 

(p-l)p^  -  0  ,  (65a) 

or 

(p-l)p‘^  -  0  ,  (65b) 

when  (h/c)  =  0  .  The  root  p  =  1  of  equations  (65)  corresponds  to 
the  exact  solution,  whereas  the  remaining  roots  p  =  0  are  the  spurious 
solutions.  For  (h/e)  =  0  ,  the  spurious  solutions  for  any  Adams  formula 
do  not,  in  fact,  exist.  Thus,  the  Adams  formulas  are  optimum  in  terms 
of  damping  out  the  spurious  solutions  when  (h/e)  -  0  .  Because  of  this 
Important  property,  the  Adams  method  is  probably  the  most  widely  used 
of  all  predictor-corrector  procedures.  By  comparing  equation  (61)  with 
(65a)  >  the  following  values  for  the  constants  A^^  are  obtained: 

A^  -  1  ,  A^  -  Ag  =  ...  =  Aj^,  =  0  .  (66a) 


By  comparison  of  equations  (62)  and  (63)  with  (65b),  and  equation  (64) 
with  (65a) ,  the  following  values  for  the  constants  a .  are  obtained 
for  all  three  cases: 


a  *  1  ,  a.  *  a^  —  ...  —  a.  —  0  . 
0^12  k 


(66b) 
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The  preceding  result  means  that  both  of  the  summations  farthest  to  the 
left  In  formulas  (48)  and  (49)  can  nov  be  replaced  by  /(xg)  •  Since 
conditions  (66b)  ailso  satisfy  the  first  of  eqxiatlons  (60),  the  latter 
Is  no  longer  needed. 

To  determine  uniquely  the  remaining  k'-f  1  constants  B.  In 
formula  (48)  by  the  p  linear  equations  (60),  k'  must  equal  p  -1  . 
It  Is,  of  course,  possible  for  k*  to  be  greater  than  p-  1  ,  thereby 
providing  k'  -  p  1  free  paraneters,  which  might  then  be  chosen  to 
Improve  stability  when  (h/c)  >  0  .  Such  an  extension  Is,  however, 
beyond  the  scope  of  this  work,  and  Is  not  given.  To  determine  uniquely 
the  remaining  k  2  constcuits  b^  In  formula  (49)  by  the  p  linear 
equations  (60),  k  must  equsd  p-2  .  The  same  remark  as  above  con¬ 
cerning  free -parameters  also  applies  to  the  corrector  (49).  The 
results  of  these  ccdculatlons  are  summarized  for  p  >  2,  3,  4  ,  euid  5 
in  the  following  table: 


CONSTANTS  B,  AND  b|  FOR  THE  ADAMS  METHOD 


p 

®o 

®1 

®2 

®3 

®4 

B 

B 

B 

^2 

^3 

3 

2 

1 

■  2 

- 

- 

1 

2 

1 

2 

- 

- 

- 

3 

23 

12 

16 

“  12 

- 

• 

8 

12 

1 

"  12 

- 

- 

4 

s 

-i 

i 

- 

i 

‘  It 

1 

- 

5 

1901 

720 

2774 

720 

2616 

720 

1274 

720 

720 

2^1 

720 

646 

720 

264 
"  720 

106 

720 

.22. 

720 

Since  k'  ■  p  -  1  and  k  ■  p  -  2  ,  the  three  relevant  character¬ 
istic  polynomleds  (6l),  (62),  and  (64)  now  tedce  the  form 


,(pr) 


P  P"1 
P  -  P 


0  , 


(67) 


1-0 
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,(co) 


.  pi>-®  ^  I  ^  bjpP-2-1 .  0  , 
1-0 


(68) 


.  pp+  (-i+b_3^  ly  b^pP"^“^- 


eZ_, 

1-0 


Z_, 

1-0 


vhere  and  bj^  eure  given  by  the  preceding  table.  For  a  given  value 
of  p  ,  each  root  of  each  of  the  characteristic  polynomials  (6t)>  (68), 
and  (69)  depends  on  the  value  of  the  parameter  h/e  .  For  a  sufficiently 
small  positive  value  of  h/e  all  of  the  roots  are  located  In  the  Inte¬ 
rior  of  the  unit  circle  In  the  complex  plane.  In  this  clrcximstance,  as 
mentioned  earlier,  the  particular  formula  Is  stable.  Thus,  the  predictor 
Is  referred  to  as  stable  when  all  of  the  roots  of  equation  (67)  sure  In 
the  Interior  of  the  unit  circle.  The  Adams  method  Itself,  Is  stable, 
of  course,  only  when  the  roots  of  equation  (69),  the  characteristic  equa¬ 
tion  for  the  overall  procedure,  are  In  the  Interior  of  the  unit  circle. 

For  Boae  positive  value  of  h/e,  equation  (6T)  will  have  a  root 
with  unit  magnitude,  l.e.,  |p|  =  1  .  When  this  Is  the  smallest  positive 

value  of  h/e  for  which  this  occurs,  then  this  value  Is  denoted  by  the 
stability  parameter  b(p^)(p)  .  Hence,  the .predictor  Is  stable  for  all 
values  of  h/e  that  satisfy  0  <  (h/e)  <  B'^^'  .  In  a  similar  fashion, 
the  stability  parameters  B'®°)(p)  and  b(P’^)(p)  are  defined. 

By  means  of  Wllf 's  (1959)  criterion,  the.  stability  parameters  were 
found  for  veurlous  values  of  p  .  The  following  table  summarizes  the 
results  of  these  lengthy  calculations: 


STABILITY  PARAMETERS  RELEVANT  TO  THE  ADAMS  METHOD 


p 

B<‘''>(P) 

2 

1.00 

00 

2.00 

3 

.545 

6.00 

1.72 

u 

.300 

3.00 

1.28 

5 

.163 

1.84 

.947 
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The  preceding  table  shows  that  all  three  stability  paraneters 
decrease  with  Increasing  p  •  This  result  Is  directly  opposite  to  that 
found  for  the  Runge-Kutta  procedure.  As  mljdit  be  anticipated,  ^9^  ? 
given  value  of  p  ,  the  overall -procedure  stability  parameter 
has  a  value  between  and  Bv®®'  .  Nevertheless,  It  differs 

appreciably  from  b(®°)  .  Thus,  Chase's  assertion  is  valid  here.  Since 
BVPc)  does  not  vary  much,  there  Is  little  reason  to  prefer  the  low- 
order  (  p  ■  2  or  3  )  Adams  procedures  to  the  hlj^r-order  ones. 

The  foregoing  analysis  Indicates  that  stability  Is  the  central 
problem  when  stiff  equations  are  Integrated  by  a  predictor -corrector 
procedure.  For  a  given  value  of  p  ,  however,  the  Adams  method  does 
not  have  em  optimum,  l.e.,  a  maximum,  stability  parameter.  We  shall 
now  consider  briefly  the  question  of  a  procedure  with  an  optimum  sta¬ 
bility  parameter.  The  discussion  is,  however,  restricted  to  p  -  4  . 

Set  k'  »  2  and  k  -  2  in  the  Mnerallzed  predictor-corrector 
procedure  given  by  formulas  (48)  and  (49) •  These  formulas  then  become 


z(x^+  h)  -  A^y(x^)+  A^y(x^-  h)+  A2y(x^-  2h) 

+  h[B^y'(x^)+  B^y’(x^-  h)+  B2y'(x^-  2h)j  , 


(70) 


y(x^+  h)  =  a^y(x^)+  a^y(x^-  h)+  a2y(x^-  2h) 

+  h^b_^z'(x^+  h)+  b^y’(x^)+  b_^y’(x^-  h)+ 

where,  as  a  result  of  equations  (60),  the  various  constants  satisfy  the 
relations 


1  -  aj^-  ag 

A^  »  -  8  -  Ag  , 

^(9  -  a^) 

»  9  , 

^(19  +13a^+  Sag) 

^(-  5  +  13a^+  32a2)  , 

■  Ic? 

^(l-ai+  8^2^  » 

®2  *  3^" 

\  (72) 
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Formulas  (70)  and  (71)  Involve  three  free  parameters  ai ,  »  and  > 

which  will  be  chosen  later.  The  corrector,  formula  (7lT>  been 
examined  in  detail  for  relative  stability  by  Hamming  (19^).  Crane  and 
Lambert  (1962)  have  also  examined  this  corrector  for  stability.  They 
find,  for  exaiiq;)le,  that  formula  (71)  is  stable  for  0  <  (h/c)  <  lOU 
when  a^  -  -  1.6  and  ■  *925  . 

Stability  of  the  overall  procedure  is  determined  by  the  roots  of 
the  characteristic  polynomial 


-  [“2-  <>>2*  ‘'.1*2>(!)*  ■  0  '  (T3) 


which  is  readily  obtained  from  equation  (614-).  Each  set  of  values  for 
ai,  &2  >  ^  results  in  a  specific  pol^omial  (73)  and  thus  a 

specific  stability  peurameter.  It  has  not  been  possible,  however,  to 
determine  analytically  the  value  of  a^,  Sg  ,  and  Ag  such  that  the 
overall  stability  parameter  Is  an  optimum.  Consequently,  stability 
parameters  for  a  wide  choice  of  values  for  a^,  ag  ,  and  were 
computed.  From  among  the  different  values  chosen,  a,  =  .8,  a2  =  0  , 
and  A2  =  0  resulted  in  the  largest  value  of  the  stability  parameter. 
With  this  choice,  the  overall  procedure  is  stable  for  0  <  (h/e)  <  1.82  . 
While  this  choice  probably  does  not  correspond  to  the  optimum  one,  the 
author  believes  that  it  does  not  differ  from  it  appreciably.  Compared 
to  the  p  z  4  Adams  method,  this  choice  results  in  a  4o^  Increase  in 
the  value  of  the  stability  parameter.  Consequently,  this  procedure 
rather  than  any  of  the  Adams  proced\ires,  will  be  used  in  the  compara¬ 
tive  evaluation  given  in  Chapter  7*0. 

6.3  CURTIS$.HIRSCHFELDER  PROCEDURE 

This  chapter  concludes  with  an  analysis  of  the  Curtiss -Hlrschf elder 
procedure.  Their  procedure  utilizes  the  following  first-order  corrector: 

y(x^+  h)  -  y(x^)  +  hy'(x^+  h)  .  (74) 


They  avoid  the  use  of  a  predictor  by  noting  that  G  and  €  in  equation 
(2a)  are  functions  only  of  x  ,  rather  than  of  x  and  y  .  In  this 
case,  erfter  equation  (2a)  is  used  to  replace  the  unknown  derivative 
y'(x  +  h)  ,  equation  (74)  can  be  solved  explicitly  for  y(x  +  h)  . 

Whili  this  is  possible  with  eqiiatlon  (2a),  it  is  generally  not  possible 
with  the  chemical  or  vibrational  rate  equations. 
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Curtlas  and  Hlrachfelder  are  concerned  with  both  positive  and 
negative  values  for  h/e  In  their  article,  even  though  the  exact 
solution  diverges  from  the  asymptote  (see  Sketch  2)  when  the  value  Is 
negative.  They  assert  that  the  numerical  solution  of  equation  (2a) 
as  obtained  by  their  procedure  will  always  converge  towcurd  the  asymp¬ 
tote  for  both  positive  and  negative  values  of  h/c  .  This  assertion 
Is  based  on  an  intuitive,  but  Incorrect,  geometric  argument.  They 
then  verify  this  assertion  by  Integpratlng  equation  (2a)  first  with 
(h/c)  ■  5  and  then  with  (h/c)  >  -  5  .  In  both  cases,  the  numerical 
solution  does  Indeed  converge  to  the  asymptote.  To  further  check 
their  assertion,  we  now  examine  the  stability  of  formula  (7^)  by 
means  of  Its  characteristic  polynomial,  which  Is  given  by 

Equation  (73)  represents.  In  a  p,(h/e)  plane,  a  hyperbola  with  one 
branch  crossing  p  >  1  and  the  other  crossing  p  *  -  1  .  Since 
|p|  >  1  for  >2<h/e<0  ,  formula  (7^)  Is  unstable  In  this  region. 

This  Instability  Is  readily  verified  by  attenptlng  to  Integrate  equa¬ 
tion  (2a)  with  a  value  of  h/c  in  this  range.  The  numerical  solution 
Is  found  to  oscillate  wildly  In  this  case,  thereby  contradicting  the 
Curtiss -Hlrschf elder  assertion.  For  h/e  >  0  ,  the  procedure  Is  always 
stable  and  the  numerical  solution  does  converge  to  the  exact  solution, 
which,  eifter  the  transient  region.  Is  essentially  the  asymptote.  When 
h/e  <  -  2  ,  the  characteristic  root  is  given  by  the  brisuich  of  the 
hyperbola  crossing  p  »  -  1  .  In  this  situation,  the  numerical  solu¬ 
tion  converges  toward  a  spurious  solution.  This  spurious  solution 
happens  to  be  the  asymptote;  the  exact  solution  actually  diverges  expo¬ 
nentially  fast  from  the  asymptote,  as  noted  In  Chapter  U.O.  The 
reason  that  the  usual  stability  condition  Is  not  applicable  here  Is 
that  the  characteristic  polynomial  (73)  has  only  one  root  but  two  branches. 
Of  the  two  breaches,  only  the  one  crossing  p  =  1  corresponds  to  the 
exact  solution. 


7.0  CONCLUDING  REMARKS 

7.1  COMPARATIVE  EVALUATION  OP  THE  DIFFERENT  INTEGRATION  PROCEDURES 

Many  numerical  procedures  are  available  for  the  Integration  of 
systems  of  ordlnsuy  differential  equations.  Basically,  three  factors 
should  be  considered  as  criteria  for  deciding  which  method  Is  best 
suited  to  a  peurtlcular  problem  -  namely,  truncation  error,  stability, 
6md  computation  time.  The  iMt  factor  Is  of  special  Importance  wlien 
the  Integration  step  size  Is  very  small  In  compco'lson  with  the  total 
distance  to  be  covered.  With  reference  to  these  three  factors,  two 
different  Integration  procedures  can  be  compeured  as  follows:  For  a 
given  step  size,  the  hlgher-ordei-  procedure  Is  more  precise.  When  two 
different  procedures  are  of  the  same  order  their  truncation  errors  are 
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here  considered  to  be  approximately  equal.  Stability  is  compared  on 
the  basis  of  B  ,  except  that  (3/U)B  is  used  for  a  predictor-corrector 
procedure  (see  equation  (46)).  Computation  speed  is  conqAred  by  the 
number  of  derivative  evaluations  per  integration  step  and  the  maximum 
step  size  as  dictated  by  stability  requirements. 

As  noted  in  Chapter  5«0,  the  fourth-order  procedure  appears  to  be 
the  most  suitable  of  the  Runge-Kutta  procedures  for  the  solution  of  the 
nonequilibrium  equations.  This  conclusion  is  based  on  the  observation 
that  the  third-order  Runge-Kutta  procedure  is  not  as  stable  or  as 
accurate  as  the  fourth-order  procedure^  whereas  the  fifth-order  pro¬ 
cedure  requires  too  much  computation  time. 

All  of  the  predictor-corrector  procedures  discussed  in  Chapter  6.0 
require  the  evaluation  of  two  derivatives  per  integration  step.  Thus, 
these  methods  differ  primarily  in  the  order  of  their  truncation  error 
and  in  their  stability.  When  these  two  factors  are  considered,  the 
method  given  at  the  conclusion  of  Section  6.2  appears  to  be  somewhat 
superior  to  the  Adams  method. 

A  comparison  of  the  fourth-order  Runge-Kutta  procedure  with  the 
fourth-order  predictor-corrector  procedure  given  at  the  end  of  Section 
6.2,  reveals  that  the  Runge-Kutta  procedure  is  more  stable  (  B  =  2.8 
vs.  (3/U)B  »  1.36  ),  but  that  the  predictor-corrector  procedure  is 
slightly  faster.  The  latter  requires  two  evaluations  of  the  deriva¬ 
tives  per  integration  step,  whereas  the  fourth-order  Runge-Kutta  pro¬ 
cedure  requires  six.  On  an  overall  basis  then,  the  fourth-order  pre¬ 
dictor-corrector  procedure  appears  to  be  slightly  superior  to  the 
Runge-Kutta  procedure. 

Despite  the  foregoing  theoretical  advantage,  the  Runge-Kutta 
procedure  was  chosen  for  the  numerical  calculations  given  in  Part  II 
because  its  simplicity  makes  it  easy  to  program  for  a  digital  computer. 
Had  the  computation  time  using  the  Runge-Kutta  procedure  proven  to  be 
unreasonably  long,  then  a  predictor-corrector  procedure  would  have 
been  tried.  The  already  programmed  Runge-Kutta  procedure,  however, 
would  still  be  needed  to  start  this  procedure  and  to  compute  inter¬ 
vening  points  if  the  step  size  is  halved.  Since  the  conqiutatlon  time, 
when  the  Runge-Kutta  procedure  was  used,  turned  out  to  be  acceptably 
short,  no  other  procedures  were  attempted. 

7.2  CONTROL  OF  THE  INTEGRATION  STEP  SIZE 

It  is  the  purpose  of  any  method  that  controls  the  step  size  h 
to  maximize  this  quantity  while  maintaining  stability  and  keeping  the 
truncation  errjr  small.  Since  eui  estimate  of  the  truncation  error 
for  a  given  step  Is  available,  the  error  condition  Is  readily  met. 

The  stability  condition  requires  that  h/e  be  slightly  less  than  the 
appropriate  stability  parameter.  The  value  of  the  quantity  €  , 
although  considered  a  constant  In  the  foregoing  analysis,  actually 
changes  continuously;  and  as  demonstrated  in  Chapter  3.0,  It  Is  not 
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easily  found.  It  would  therefore  be  desirable  to  develop  a  nethod 
for  controlling  h  that  does  not  require  the  evalxuitlon  of  h/e  . 

The  only  readily  available  quantity  for  controlling  h  Is  the 
estlaated  truncation  error  per  step.  For  exaniple,  Runge-Kutta  Inte¬ 
gration  generally  controls  the  step  else  by  a  technique  that  attenq>ts 
to  hold  constant  the  estlaated  error  per  step.  As  discussed  In  Part 
II,  this  method  results  In  cm  erratic  variation  of  step  size  when  the 
rate  equations  of  Part  II  are  Integrated.  In  terms  of  maximizing  the 
step  size  and  still  maintaining  stability,  this  method  Is  poor. 

A  method  that  is  based  solely  on  the  estimated  truncation  error 
for  a  single  step  is  not  feasible,  since  there  Is  no  direct  connection 
between  the  magnitude  of  this  error  and  the  stability  of  the  Inte¬ 
gration  procedure.  It  is  possible,  for  example,  for  a  stable  procedure 
to  have  a  large  error  in  a  given  step  or  for  an  unstable  procedure  to 
have  a  small  one.  The  stability  of  the  procedure,  however,  determines 
the  rate  of  growth  of  the  error  over  successive  steps.  When  a  pro¬ 
cedure  Is  unstable,  the  error  will  Increase  with  each  succeeding  step, 
whereas  the  reverse  Is  true  when  the  procedure  Is  stable.  The  Improved 
technique,  given  In  detail  in  Part  II  for  controlling  the  step  size 
when  the  Runge-Kutta  procedure  is  used  is  based  on  this  concept. 

In  this  method,  three  constants  R^,  R2  ,  and  R^  are  defined  by 

R^  ■  lower  reference  value  for  the  estimated  truncation  error, 

Rg  *  upper  reference  value  for  the  estimated  truncation  error, 

R^  a  number  of  integration  steps  between  possible  increases  in  h  . 

These  constants  approximately  measure  the  rate  of  growth  of  the  trun¬ 
cation  error.  When  the  estimated  truncation  error  is  less  than  R^ 
for  R^  consecutive  steps,  the  procedure  is  then  considered  suffi¬ 
ciently  stable  to  allow  for  a  small  Increase  in  step  size.  On  the 
other  hand,  if  the  estimated  truncation  error  grows  with  each  succeed¬ 
ing  step  such  that  it  ultimately  exceeds  R2  ,  where  0  <  R^  S  Rg  , 
then  the  procedure  is  considered  unstable  euid  the  step  size  is  reduced 
by  a  small  amount.  The  percent  Increase  or  decrease  in  step  size  is 
governed  by  two  constants.  To  avoid  going  from  a  stable  to  an  unstable 
integration  condition,  the  percent  Increase  in  h  should  be  small, 
e.g. ,  30^.  Similarly,  to  avoid  going  from  a  slightly  unstable  condition 
to  an  overly  stable  condition  that  would  require  excessive  computation 
time,  the  percent  decrease  in  h  should  also  be  small,  e.g.,  30^.  As 
verified  by  Figure  1  of  Part  II,  the  above  method  approximately  maxi¬ 
mizes  h  while  at  the  same  time  stability  is  maintained. 
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APPiNPIX  I 

PROOP  THAT  ALL  OP  THE  EIGENVALUES  <<1  ARE  REAL  AND  POSITIVE 


This  appendix  proves  the  following  fundamental  theorem  referred 
to  In  Chapter  3.0: 

Theorem:  All  eigenvalues  (i=l,...,N5)  ,  given  by  the  chemical 

characteristic  equation 

-  «ul  =  0  .  (2'*») 

are  real  and  positive. 

A  sequence  of  five  lemmas  is  used  to  establish  this  theorem.  A 
few  preliminary  remarks  concerning  the  nature  of  the  proof  will  be 
helpful. 

Only  two  physical  concepts  are  required  in  the  proof.  First,  the 
equations  that  express  conservation  of  components 


n 


i 


i  =  N3+1,...,N^  ,  (l-l) 


are  necessary.  Second,  the  structure  of  the  chemical  rate  equations (as 
given  in  Chapter  3.0),  when  the  flow  is  close  to  equillbrivm,  is  also 
involved.  Thus,  the  proof  does  not  require  conservation  of  momentum 
or  of  energy.  The  theorem  is  therefore  valid  for  Inviscid  or  viscous 
flows  as  well  as  non-radiating  or  radiating  flows.  Likewise,  the  second 
law  of  thermodynamics  is  not  required  in  the  proof.  The  theorem  as 
given  here,  is  not  applicable  when  certain  physical  processes  that 
alter  the  form  of  the  chemical  rate  equations,  such  as  diffusion,  are 
to  be  considered. 

Two  matrices,  {(kj^^))  and  ((^wf))  >  whose  elements  are 
and  (see  equations  (l6)  and  (l8)),  are  used  throughout  the  proof. 

The  first  is  an  N3  by  Nt  matrix,  while  the  second  is  an  by 

matrix.  Although  the  theorem  is  directly  concerned  only  with  the  N3 
eigenvalues  of  ((^ki))*  determined  by  equation  (2Ub),  most  of  the 
proof  actually  deals  with  ((kj^^))  .  The_basic  reason  for  this  is  the 
greater  simplicity  of  the  elements  of  ((«>/))  as  compared  with  those 
of  ((Ki,,))  . 

We  here  follow  standard  matrix  theory  terminology  by  referring  to  the 

eigenvalues  determined  by  equation  (24b)  as  belonging  to  the  matrix 
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The  five  leonas  can  be  briefly  sunsnarlzed  as  follows:  LensBa  1 
establishes  two  useful  relations  necessaxy  In  the  subsequent 
development.  Leimna  2  establishes  the  connection  between  eigenvalues 
of  ((k]^^))  and  those  of  ((kj^^))  •  Leimna  3  demonstrates  that 
((kj^^))  has  N3  non-zero  eigenvalues.  Lemmas  4  and  ^  Eu:e  used  to 
prove  that  la  positive  semldef Inlte.*  Once  these  lemmas  are 

established,  the  actual  proof  of  the  theorem  then  follows  readily. 
Following  this,  the  appendix  concludes  with  the  demonstration  of  two 
corollaries  that  help  clarify  the  preceding  proof.  The  first  of  these 
demonstrates  that  It  Is  Immaterial  which  side  of  a  reaction  Is  taken 
to  refer  to  the  reactants  and  which  to  the  products.  The  second  shows 
that  the  eigenvalues  of  ((•‘iti))  do  not  depend  on  the  specific  choice 
of  the  mole-mass  ratios  ti=N3+l,. . .  ,N-]^)  that  are  taken  to  be 

given  by  the  conservatlon-of -components  equation  (l-l).  Both  of  these 
corollaries  could  have  been  anticipated  from  purely  physical  reasoning. 


Lemma  1: 


(a)  The  quantities  are  given  by 


'll 


1  V  rrri 

L  et  * 


r-1 


1  -  1,...,N^,  i  •  1,...,Nt  .  (1-2) 


(b)  The  relations 


”3 

'u  ■  Z  “ij'ji  • 

>1 


1  .  Nj+  l,...,Nj^,  i  ■  l,...,Nj^  (1-3) 


are  a  consequence  of  the  equations  expressing  conservation  of  conqponents 

(l-l). 


A  matrix  is  referred  to  as  positive  semidefinlte  if  all  of  its  eigen¬ 
values  are  real  and  non-negative.  If  all  the  eigenvalues  of  a  matrix 
are  real  and  positive,  l.e.,  none  are  zero  valued,  the  matrix  is 
called  positive  definite. 
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Proof:  Differentiating  equations  (6b)  gives 


3_-a 

Mij: 

.(T) 


cr' 


P 


-  (‘■r-  1)  ^ 
r  n, 


r  ■  i  ■  f  (l“U) 


where  p,T  ,  and  all  nj^  ,  except  n^  ,  are  held  fixed  in  accordance 
with  equations  (ll).  Since  (Ly)*  ■  0  (see  equation  (8a)),  the  refer¬ 
ence  value  for  SLy/dnj  is  therefore  given  by 


r  *  1,«<.,N2  t  ■  1, . . . ,N^  .  (l"5) 


Substituting  equation  (l-5)  into  equations  (l6)  results  in  equations 
(l-2),  thereby  establishing  part  (a)  of  the  lemma.  It  is  important  to 
note  that  equations  (l-5)  are  valid  only  if  the  reference  state  is  given 
by  the  equilibrium-flow  solution  or  the  local  equilibrium  value. 

Part  (b)  is  proved  by  differentiating  equation  (l-l)  and  then  re¬ 
placing  Dnj/Dt  by  equations  (l5b)  as  follows: 


i 


I 

! 


J 


r  (n 


l,e 


np  , 


i  -  N3+1,...,N^  .  (1-6) 
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Tba  i  and  J  avinmatlons  have  been  Interchanged  in  the  last  step. 

"By  comparing  this  result  vith  equations  (I5b),  we  obtain  equations  (l-3). 
Thus,  the  lenma  Is  proved. 

Leosaa  2  vlll_prove  that  any  eigenvalue  of  ((K]ci))  Is  also  an 
eigenvalue  of  ((kj^j))  •  In  addition,  the  lenma  will  show  that  any 
eigenvalue  of  ((kj^i))  that  Is  not  an  eigenvalue  of  ,  Is  zero. 

These  two  results  are  established  by^dewnstratlng  that  ((kiij))  and 

((kj^j))  ,  except  for  the  factor  (a)  ^  ^  ,  have  Identical  characteristic 
eqviatlons. 

Lenma  2;  The  relation 

N  •N  *  * 


Is  valid. 


Proof;  The  left-hand  side  of  equation  (l-7)  is  a  polynomial  In  k  of 
degree  •  This  polynomial  Is  denoted  by  P  ■  P(k)  ,  and  when  expli¬ 
citly  written  out  is  given  by 


K 


'll 


‘12 


'13 


P(k) 


K  -  K 


?2 


-  a, 


23 


•  (l-&a) 


R  -  K 


NiNi 


First,  column  N3+I  of  P  Is  multiplied  by  and  the  result 

Is  added  to  column  1.  Commas  are  now  used  between  certain  double  sub¬ 
scripts  such  as  N3+I  and  1  ,  solely  for  clarity.  Next,  column 
N3  +  2  is  multiplied  by  ajj  ^  and  again  the  result  Is  added  to 

column  1.  By  repeating  this  process,  we  finally  obtain  for  the  first 
column  of  P 
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K  -  K 


11 


'21 


-  K. 


V 


■  ■‘n.+I,!’^  “  ■  Z  '^N.+l,; 

"  ®  J.H3+I  ’ 

»1 

■  \+2,l  ■  Z  X+2,^ 

^  ^  J-N3+I  * 

•V,i*\a'‘’  L 

^  >N3+1  ^ 


\  (1-9) 


where  equations  (I8)  are  used  to  obtain  the  top  N3  elements.  Column 
2  Is  next  altered  In  a  slmlleur  manner,  l.e.,  by  multiplying  column 

N3+lof  P  by  aji  2  adding  the  result  to  column  2  of  P  ,  etc. 
3 

After  column  2  the  process  Is  repeated  for  columns  3  through  of 
P  ,  thereby  resulting  In 
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ColuonB  N3+I  throu^  %  of  (l-8b),  of  course,  are  Identic^  with 
the  corresponding  columns  of  (l-8a).  The  top  N3  rows  of  (I-ob)  are 
now  multiplied  by  -  aj^  (lo«l,. . .  ,N3)  and  the  result  added  to 

the  H3 +1  row  of  (l-8b).  The  first  element  In  this  row  then  becomes 


>H,+1 


‘j,l 


'll 


♦  •113+1,2  •‘21 


+  . . .  + 


•li3+l,N3  “Mj,! 


■  -  K, 


H3+I 


,1 '  L  *j,i  V^' 

J-H3+I 


itwi 

h  kl 


(l-lOa) 


When  is  replaced  by  equations  (18),  and  the  result  Is  regrouped, 

we  obtain 


According  to  equations  (l-3),  each  term  enclosed  by  parentheses  In 

(l-lOb)  Is  zero.  Consequently,  the  first  element  of  row  N3+I  of . 

(l-8b)  is  zero.  In  a  similar  fashion,  each  element  from  the  first  to 
the  of  this  row  of  (l-8b)  Is  shown  to  be  zero.  The  N3+I 

element  of  this  row  becomes 


K 


+ 


*11,1*3+1 
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Which,  by  virtue  of  equations  (l>3),  is  equal  to  k  .  In  a  similar 
fashion,  the  remaining  elements  of  the  1  row  of  (l-8b)  are  shown 
to  be  zero.  Hence,  this  row  contains  only  zero  elements,  except  for 
the  diagonal  element,  which  Is  k  .  The  argument  for  the  other  rows 
of  (l-8b)  from  row  N3  +  2  through  row  is  the  same.  Thus,  deter¬ 
minant  (l-8b)  can  be  written  as  follows: 


1 


I 


t 


P  - 


s  -  It 


11 


-  K 


21 


-  K 


N-,1 


-  R 


12 


-  R 


R  -  R, 


22 


•  '‘H3,2  •••'*■ 


1,1*3  * 


.  .  .  *  H 


1,H, 


-  R. 


"N3,Hi 


K  0  ...  0 

0  K 

» 

•  •  • 


0  »  s  *  0  K 

(l-8c) 


This  determinant,  however,  is  precisely  equal  to 


N,  -N_ 


which  proves  the  lemma. 

The  purpose  of  the  next  lemma  is  to  show  that  ((it]^/))  has  N3 
non-zero  eigenvalues.  Since  ((kj^i))  Is  an  N3  by  N3  matrix,  it 
is  sufficient  to  prove  that  the  reuik  of  is  N3.*  Implicit 

in  the  proof  of  this  lemma  is  the  assumption  that  No  2  N3  .  As  noted 
in  Part  II,  an  immediate  consequence  of  this  assumption  is  that  the 
minimum  nu^er  of  chemical  rate  equations  is  N3  . 


If  the  rank  of  ((kj^/))  is  N3  ,  then  the  determinant  of  ((thi)) 
is  not  zero  and  the  matrix  hsis  no  zero  eigenvalues.  If,  however,  the 
rank  of  ((kic/))  is  less  than  N3  ,  then  the  determinant  of 
is  zero  and  the  matrix  has  at  least  one  zero  eigenvalue. 
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Lenaaa  3;  The  rank  of  ((kj^j))  Is  Nj  . 

Proof!  If  the  rank  of  ( (kj^j) )  were  less  than  N3  ,  say  N3-  1  ,  then 
there  would  exist  a  non-zero  constant  vector  (c^^)  (iBl,...,N3)  such 
that 


"3 

) 

1-1 


I  »1«U  -  0  , 


JL  *  S  •  • 


(l-ll) 


Equations  (17)  could  thus  he  transformed  as  follows: 


»3  »5  K 

I  '1  dt  •  I  ”j)  L  w,  ■ 

1-1  Ul  1-1 


{l-12a) 


The 


o 

terms  ^  equal  to  zero  as  a  result  of  equations  (l-ll). 

1«1 


Equation  (l-12a)  could  therefore  be  Integrated  to  yield 


I  'i"i 
1-1 


constant. 


(I-I2b) 


This  result  Is  not  possible,  however,  since  equations  (17)  represent 
the  minimum  possible  number  of  rate  equations.  Consequently,  the  rank 
of  inust  be  N3  and  the  lemma  is  established. 

The  remaining  two  lemmas  are  used  to  show  that  ((K]^i))  Is  positive 
semldeflnlte.  A  sequence  of  determinants,  given  by 


•'ll  **12  "'im 


« 


m  —  1 , . • « , 


(1-13) 
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forms  the  starting  point  for  the  proof.  Basic  to  the  proof  is  the 
following  theorem  from  matrix  theory;  An  by  matrix  ((kj^i)) 
is  positive  semldefinlte  If^  and  only  if, 

a  0  ,  m-l,...,Nj^  .  (I-14) 


The  above  theorem  is  usually  proved  for  positive  definite  matrices.  In 
this  case,  the  symbol  (  a  )  in  conditions  (l-*l4)  should  be  changed  to 
the  symbol  (  >  )•  A  proof  of  this  theorem  can  be  found  in  Frazer,  Duncan, 
and  Collar  (1947). 

Our  ultimate  purpose  is  to  show  that  TOsltive  definite. 

This  goal  will  be  achieved,  in  part,  by  proving  positive  seml- 

deflnlte.  This  round-about  approach  is  undertaken  because  a  direct 
proof  that  ((kj^^))  Is  positive  definite  would  be  considerably  more 
difficult  owing  to  the  complicated  structure  of  the  elements  of  ((kjj/)). 

The  actual  structure  of  the  chemical  rate  equations,  when  the  flow 
is  close  to  equilibrium,  is  embodied  in  the  quantities  ,  as  given 

by  equations  (1-2).  This  structxire  enters  in  an  essential  way  into  the 
next  two  lemmas.  To  put  it  differently',  ((i^ici))  positive  semi- 
definite  because  of  this  structure. 

Before  starting  the  next  lemma,  we  will  define _certain  quantities 
that  are  necessary  for  its  formulation.  When  the  ,  given  by 

equation  (l-2),  are  substituted  into  equation  (l-13),  the  elements  of 
each  column  are  seen  to  contain  the  factor  l/n^  .  Consequently,  a  new 
sequence  of  determinants  is  defined  as  follows:  ■ 


(l-15a) 


I 

1^1 


''rl^l 

r 


r"l 


V  V 

rm  rl 
r 


I 

r«l 


'^rl'^rm 

r 


r«l 


V  V 
rm  rm 
-gs— 
r 


,  iii^l,...,N^  .  (l-15h) 
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When  the  various  elenents  of  determinant  are  multiplied  against 

each  other.  In  accord  with  the  usual  rules  for  congutlng  the  vadue  of 
a  determinant,  certain  quantities,  such  as  (in  this  Illustration  m  -2  ) 


(eif 


1  1  1 


will  reguleurly  appear. _.TteBe  quantities  are  here  considered  to  be  the 
variables  upon  which  depends.  The  coefficients  of  these  variables 

are  readily  seen  to  be  functions  only  of  the  consteuits  .  In  order 

to  introduce  a  notation  for  these  variables,  a  subscripted  r,  rg  ,  is 
now  used  to  denote  any  value  of  r  from  1  to  N2  •  For  a  fixed  m  , 
the  above  variables  cure  defined  by 


H(m;rj^,r2,...,r  )  s  g  } — ^  ,  (I-I6) 

Rj 


where  the  dependence  of  H  on  the  parameters  m  and  r^,...,rj^  is 
explicitly  Indicated.  For  exaii5)le,  with  m  *  2,  r^^  =  5  ,  and  r2  *  1  , 
H  Is  given  by 


H(2j5,1)  -  H(2;1,5) 


(1-17) 


For  reasons  that  will  become  apparent  in  the  proof  of  Lemma  4,  the  defi¬ 
nition  of  H  given  above  is  restricted  to  prevent  any  Q*  from  appear¬ 
ing  more  than  once.  Thus,  in  addition  to  definition  (l-lB),  we  also 
require  that  if  s  ^  t  then 


r 

s 


Different  H's  will  be  distinguished  by  subscripts.  Therefore, 
H2^(m;r|^^ . ,r^^^ )  must  have  at  least  one  r^^'  value  that  is  differ¬ 
ent  from  all  the  r(^)  values  in  H2(m;r^^' , . . . ,r^^) )  . 

Lemma  4:  Each  m  by  m  determinant  can  be  written  as 
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■  Z  . r^‘>)  , 

i-1 


m 


1,...,N^  .  (1-18) 


where 


N^l 


(m^)  (Ng-  m)l  ml  ' 


(1-19) 


and  the  constants  =  Aj^(ni;r|^\. . .  depend  only  on  the  para¬ 

meters  Vpj  . 

A  simple  Illustration  of  Lemma  U  Is  provided  when  m  =  1  .  In 
this  case,  we  have 


Ng  2 

=  n*K^°^  =  V  ^  , 

"l*^l  /_,  e*  > 

1  r 
r*l 


(1-20) 


where 


(l^)  =  N2,H^(l;rJ^h  =  (0^)’^  ,  and  A^  =  . 


i?(l) 


Proof;  Each  element  of  '  contains  N2  terms,  according  to  deter¬ 
minant  (l-15b).  Therefore,  using  the  determinant  rule 


*1 


Cl  +  Cg 


Cl  d 


^(1) 


we  can  write  '  as  the  sum  of  (Ng)  determinants.  Each  of  these 
Is  an  m  hy  m  determinant,  whose  elements  have  the  form 


^rl^rj 

r 


(1-21) 


Since  every  element  In  a  particular  column  contains  Vj,j/0*  ,  this 
factor  may  be  taken  outside  of  the  determinant.  Thus,  each  of  the 
(Ng)™  determinants  Is  of  the  form 
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'’r  «  •  *  • 


(1-22) 


Note  that  If  =»  s  ^  t  ,  the  determinant  In  expression  (l-22)  Is 
zero  since  two  columns  are  equal.  Hence,  only  those  determinants  for 
which  rg  4  s  t  need  be  considered.  In  addition,  If 

m  >  Ng  then  at  least  two  columns  must  be  Identical  In  each  of  the 
(N2)™  determinants.  Consequently,  we  have 


(1-23) 


for  all  m  >  . 

According  to  the  foregoing,  1^^)  Is  a  linear  function  of  the 
variables 


(I-2I;) 


where  rg  ’f  r^  when  s  t  .  Each  of  the  products  In  expression  (l-2U) 
therefore  contains  m  distinct  selected  from  the  No  distinct 
^  .  Thus  the  number  of  dlstljjct  variables  of  the  form  (l-2lt)  Is  given 

by  the  binomial  coefficient  (m^)  .  Each  of  these  was  defined  earlier 
as  an  H^(m;  r|^^),...,r4^')  .  Therefore,  K^^)  can  be  written  as 
equation  (I-I8),  thereby  proving  the  lemma. 
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\ 


Lemia  5t  The  constants  Aj^  are  given  by 


V. 


(v) 


^2 


1  - 


(1-25) 


Proof;  The  number  of  permutations  of  the  m  distinct  numbers 

is  ml  .  Consequently,  the  variable  will  occur  ml 

times  in  the  original  expansion  of  into  (N2)™  determinants. 

For  example,  with  N2  =  3  and  m  =  2  ,  each  of  the  (N2)®  =  (3)^  =  9 

determinants  In  the  expansion  Is  multiplied  by  one  of  the  follow¬ 

ing  variables: 


■ 


As  we  have  already  seen,  the  coefficients  of  the  (6^)"^  variables  are 
zero.  There  are  then  ml  =  21  =  2  variables,  such  as  and 

,  each  of  which  is  designated  by  the  same  . 

According  to  expression  (l-22),  each  of  the  ml  terms  in  the  coef¬ 
ficient  of  a  particular  Hj  is  of  the  general  form 
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With  fixed,  we  now  examine  the  ml  determinants  appearing  in  the 
coefficient  of  this  .  Bach  of  these  determinants  contains  the  seuae 

columns  but  in  a  diffennt  sequence.  In  fact,  no  two  of  the  ml 
determinants  have  their  columns  in  precisely  the  same  sequence.  Since 
interchanging  columns  in  a  determinant  slaqply  changes  its  sign,  all  ml 
determinants  are  equal  to 


(sffi  o) 


(1-27) 


where  sgn  a  is  1  or  -  1  depending  on  whether  the  permutation  o 
is  even  or  odd.  One  column  Interchange  represents  an  odd  permutation, 
two  an  even  permutation,  etc. 

The  constant  is  then  obtained  by  summing  the  ml  terns.ln 
the  coefficient  of  as  follows: 


where  the  summation  is  over  all  ml  permutations  of  r^^^ ,. . . ,r^^^  . 
After  the  determinant,  which  is  a  common  factor,  is  taken  outside  the 
summation,  A^  becomes 
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(l-28b) 


According  to  the  usual  definition  of  a  determinant  (see  Marcus  (i960)), 
the  summation  is  given  by 


(1-29) 


Consequently,  Is  given  by  equation  (l-25),  thus  proving  the  lemma. 

A  simple  illustration  of  Lenma  ^  is  provided  by  the  earlier  example 
when  N2  =  3  and  m  =  2  .  The  six  non-vanishing  determinants  In  the 
expansion  of  are  then  treated  as  follows: 


AE0C<TDR.«3.I2 


2 

2 

2 

"ll  ''2l''22 

V  ''ll''l2 

'’ll  ''3l''32 

W  ""g» 

^  TT 

T  "T" 

4- 

+ 

2 

.  .  2 

2 

''ll'’l2  ''22 

'^2l'^22  ''12 

'^ll''l2  ''32 

8^ 

®1  ®3 

f(i). 


^11^22 


'’ll  ''21 

'’21  '’11 

.  '’2l''l2 

'’12  ''22 

'’22  ''12 

-  ''2l''l2>‘ 

.  ^''ll'’32 

'’31 


''l2  '’32 


+  •  •  • 


12 


m 


(1-30) 


A  second  illustration  of  the  lemma  is  provided  by  ,  given  by 

equation  (l-20). 


As  a  result  of  Lemmas  4  and  the  proof  that  ((k]^^))  is  positive 
semidefinite  is  simple.  Since  all  and  A^^  are  non-negative,  all 

are  non -negative.  Equations  (l-15a)  then  imply _that  all  are 

non-negative.  Conditions  (l-l4)  then  state  that  ((k]^^))  is  positive 
semidefinite,  i.e.,  all  its  eigenvalues  are  real  and  non-negative. 


The  proof  of  the  theorem  stated  at  the  outset  of  this  Appendix  is 
now  also  simple.  The  preceding  peuragraph  demonstrated  that  all  eigen¬ 
values  of  ((K]£i))  real  and  non-negative.  According  to  Lenmia  2, 
however,  all  eigenvalues  of  ((kv/))  are  eigenvalues  of  ((itki))  • 

Thus,  all  eigenvalues  of  are  real  and  non-negative.  Lemma  3, 

however,  proved  that  all  N3  eigenvalues  of  ((kjj|))  are  non-zero. 
Consequently,  all  eigenvalues  of  ((ki^/))  are  real  and  positive,  thereby 
proving  the  theorem. 

Two  corollaries  of  the  theorem  are  now  given.  Corollsuy  1  is 
primarily  a  consequence  of  Lemma  1. 
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Corollary  1;  The  specific  values  for  the  Nj  eigenvalues  Kj^  are 
unchanged  when  the  forward  and  bachward  directions  of  any  reaction  are 
Interchanged. 

Proof;  All  that  is  necessary  is  to  show  that  the  are  unaffected 

when  the  two  sides  of  a  reaction  are  interchanged.  For  convenience, 
asBvuae  that  reaction  1  is  so  altered,  while  the  remaining  reactions  are 
left  unchanged.  All  Vi4  must  then  be  replficed  by  -  .  Also, 

according  to  equations  toe). 


(1-31) 


Noting  that  Kc;,(T»)  =  and  Li(p*,  T*,  n^)  =  0  , 

(see  equations  t^))>  we  can  write  equations  (6by  for  r  =  1  as  follows: 


1 


Pl-Ol 

WT~ 


(I -32a) 


Multiply  both  sides  of  equation  (l-32a)  by 


NU 

n 


te=N^+l 


(n*) 


'Ik 


=  1 


(1-33) 


since  Viit  =  0  for  k  =  N]^ +1,. . .  ,Ni|^  ,  and  then  rearrsmge  the  results 
to  obtain 


(p*) 


“l 


-1 


"fl 


(T*)  n  ("v) 

k^l 


a. 


Ik 


Pl-i 


=  (p*)  ^  n 


'Ik 


(l-32b) 


k=l 


By  comparing  equations  ( I -31)  and  (l-32b),  we  observe  that  (0^)  is 
unaltered  when  the  two  sides  of  reaction  1  are  interchanged.  When  these 
changes  are  introduced  into  the  r  =  1  term  of  kw  >  as  given  by 
equations  (l-2),  we  see  that  remains  unchanged^ 


Corollary  2;  The  specific  values  for  the  Nj  eigenvalues  it^  do  not 
depend  on  the  particular  choice  of  mole-fnass  ratios  nj  (iaN*+l,. . .  ,Nj^) 
taken  to  be  given  by  the  conservatlon-of -component  equations  (I-l). 
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As  an  llluBtration  of  Corollary  2,  consider  the  model  for  air  given 
In  Part  II,  where  the  mole-mass  ratios  represent  0,  NO,  , 
and  Og  ,  respectively.  From  the  two  equations. 


"1 

+  +  2n^  - 

constant 

"2 

+  nj  +  2nj^  = 

constant 

(I-3U) 


that  express  conservation  of  atomic  oxygen  and  atomic  nitrogen,  respec¬ 
tively,  eight  different  sets  of  equations,  each  set  containing  two 
equations  of  the  form  (l-l),  are  equally  possible.  The  corollary  then 
asserts  that  the  eigenvalues  are  the  same  no  matter  which  of  the  elgd^t 
sets  of  equations  (l-l)  Is  chosen.  While  the  eigenvalues  are  Invariant 
with  respect  to  this  choice,  the  coefficients  in  the  linear  trans¬ 

formation  (l9a)  do  depend  on  the  choice. 

Proof:  Corollary  2  Is  established  by  noting  that  all  eigenvalues  of 
are  also  eigenvalues  of  ((^k^))  >  as  proven  by  Lemma  2.  The 
eigenvalues  of  ((*ti5«))  ,  however,  do  not  depend  on  the  conservation  of 
conqponent  equations  tl"l)*  Consequently,  the  corollary  is  proven. 

The  eigenvalues  can  be  computed  either  by  solving  the  N3  by  N3 
determinant  equation  (24b),  or  by  solving  the  N]^  determinant 

equation 

=  0  •  (1-35) 


This  latter  equation  yields  the  correct  N3  non-zero  eigenvalues  plus 
Nj^-N3  zero  eigenvalues.  While  the  elements  of  equation  (l-35)  are 
considerably  simpler  than  those  of  equation  (24b) ,  the  higher  order  of 
equation  (l-35)  more  than  offsets  this  advantage.  Thus,  for  computa¬ 
tional  purposes  equation  (24b)  Is  superior.  This  is  especially  true 
when  the  eigenvalues  are  to  be  computed  with  a  digital  computer  using 
a  standeurd  eigenvalue  routine.  These  routines  sometimes  have  conver¬ 
gence  difficulties  when  there  are  multiple  zero  eigenvalues. 
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APPINOIX  II 

EXAMPLE  ILLUSTRATING  THE  TRANSFORMATION  TO  CANONICAL  FORM 
OP  THE  CHEMICAL  RATE  EQUATIONS 

By  means  of  a  simple  excunple,  this  appendix  Illustrates  the  theory 
given  In  Chapter  3.0.  The  gas  model,  Including  the  reactions,  Is  not 
presumed  to  he  physically  realistic;  It  was  chosen  purely  on  the  basis 
of  Its  simplicity.  In  addition  to  the  foregoing,  the  appendix  also 
demonstrates  the  behavior  of  the  transformation  given  In  Chapter  3.0  In  the 
frozen  and  equilibrium  limits.  This  Is  accomplished  by  applying  the 
limit  processes  to  the  Illustrative  example. 

We  assume  a  gas  composed  of  the  three  reacting  species  X,  X2  and 
X3  .  Their  mole-mass  ratios  are  designated  by  n^^,  n2  >  and  n^  , 
respectively.  Three  reactions  are  assumed  to  govern  the  chemical  kine¬ 
tics  as  follows: 


r  = 

1  : 

X-^X  X, 

•3  i 

r  = 

2  : 

Xg3?=2X  , 

r  = 

3  : 

X  +  X3^2X2  . 

The  gas  model  thus  contains  3  species,  1  component,  3  reactions,  and 
no  Inert  or  catalytic  bodies.  Hence,  the  parameters  are  given  by 

=  Ng  =  =  3  ,  N3  =  2  .  (Il-l) 

1 

i 

I  Conservation  of  component  X  Is  given  by 

n^  +  2n2  +  Sn^  *  3a3Q  =  constant  .  (ll-2a) 

By  solving  for  03  ,  we  obtain  the  following  equation  for  conservation 
of  components: 


*^3  ■  ^20*  ^  ^32^2  ' 


^5  “  ®30"  3  "1  ■  3  "^2  * 


(Il-2b) 
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Tbe  quantities  L,.  and  6~^  defined  by  equations  (6b)  emd  (6c) 
are  then 


Kcl  "3 


La  ■  1  -  n” 

^  ^c2  “2 


2 

1  "2 

T  s  1  • 

T  K  -  n^n. 

c3  13 


-  lSj.j_n3  , 


>  ®2  “  ^f2"2  * 


®3^  “  P^f3V3  ' 


/  (II-3) 


and  the  rate  equations  (3)  are  given  by 


^2  ^3 


^2  ^2  ^3 


(II-4) 


No  rate  equation  for  n^  is  necessary  since  this  quantity  Is  determined 
by  equation  (II -2b). 

The  coefficients  In  the  linearized  form  (l5b)  of  the  chemical 

rate  equations  are  as  follows: 
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*“  1  .  ^  1 


"2^22  "  ^  ^  ^  4  ' 


*-  1^1 
V33  *  a*  0*  " 


*-  *-  12  2 

"2'‘i2  “  V21  “  0*  ■  e*  ■  0»  ' 

X  £  o 


(II-5) 


^s'^is  “  "l'‘31  “  "  0f  0*  " 

X  0 


'^3'‘23  “  '^'*32  “  ’  ‘  4 


The  coefficients  that  accotint  for  conservation  of  con5)onents  in 


the  linearized  form 


of  the  chemical  rate  equations  are  as  follows; 


K22  =  ''22+ 


«21  -  •'ai^ 


«• 

1  , 

f  1 

1 

1  \ 

^  1  , 

/  1  1 

1  ' 

'3l'*13  " 

+ 

3 

1 

1  , 

f  1 

2 

1  \ 

^  1 

^  k 

/  1  1 

1 ' 

'32''23  “ 

0f  ' 

+ 

3 

IFJ 

0?n* 

3 

1 

'  2 

3' 

2  2 

3 

'  2 

3 

1 

(  1 

2 

1\ 

2 

2 

/I  ^  1 

1> 

■32'' 13  “ 

+ 

3 

0?n^ 

■  ^ 

te  3 

3* 

1 

'  2 

3' 

2  2 

3 

\  2 

1 

( 1 

1 

1 

2 

2  , 

1  ' 

'3l'^23  * 

Ui 

+ 

3 

■^3' 

K  '3 

ri*" 

"3‘ 

11-6) 


The  transformation  equations  (l9a)  and  (26)  are  given  by: 


qi  =  ni  +  B^2"2  > 


^  -  B^^n^  +  n^  , 
^5  ‘  ^5  ' 


(II"7) 
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where  ■  1  Is  arbitrarily  assumed.  Before  the  constants  Bj^  and 
B21  can  be  found,  the  characteristic  equation  (24b)  must  be  solved. 
This  equation  is: 


“  ■  •'ll 

**12 

“21 

K  -  K, 

t 

22 


=  0  , 


(II -0a) 


or 


K  -  Kgg)  K  +  Kg^K^)  =  0 


(II -8b) 


The  solution  for  the  two  eigenvalues  is  then  given  by 


•‘l*  2  ^11  ^*^22 


•^2=- 


■I 

i["ll 


(k11-  Kgg)  +  4k2i.c^  ,  , 


>  (II-9) 


-  +  4k2^k^2  f  , 


11  22 


where  Is  arbitrarily  chosen  to  go  with  the  solution  containing 
the  minus  sign. 

The  constants  Bjo  and  B22  are  determined  by  equations  (23). 
For  example,  Bj^  >  with  1  =  1  and  i  =  1  is  given  by 


B 


12 


or,  with  i  *  1  euid  £  *  2  ,  by 


B 


12 


(ll-lOa) 


(II -10b) 


These  two  values  for  B12  t  virtue  of  equation  (ll-8b),  are  equal. 
With  1  =  2  and  i  =  2  ,  we  obtain  B21  as 


B 


21 


*^2'  *^22 
Ki2 


(ll-ll) 
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Thus,  the  transformation  (lI-7)  is  uniquely  determined. 

The  two  transformed  rate  equations  that  replace  equations  (ZI-U) 
are  given  by  equations  (20).  They  need  not  be  repeated  here. 

Although  the  theory  of  Chapter  3.0  is  valid  only  for  a  near¬ 
equilibrium  flow^  the  frozen  limit  nevertheless  makes  physical  sense. 

For  simplicity,  ve  shall  apply  this  limiting  process  to  one  reaction; 
the  other  reactions  will,  of  course,  behave  in  a  similar  fashion  in  this 
limit.  Thus,  ve  imagine  a  situation  where  one  reaction  rapidly  changes 
from  a  nesur-equillbrlum  value  to  a  frozen  value.  To  be  more  precise, 
we  let  reaction  1  freeze  while  reactions  2  and  3  are  unaltered.  In 
mathematical  terms,  we  require  that 

-»  0  and  ,  (11-12) 


which  implies  that  (see  equations  (11-3)) 

e  and  -  Q  -  -  0  .  (11-13) 


Thus,  the  0i-term  vanishes  from  all  and  reaction  1  no  longer 

affects  thf  eigenvalues  or  the  coefficients  .  Freezing  reaction  1 

is  therefore  tantamount  to  eliminating  it  from  the  original  set  of 
reactions . 

We  next  Investigate  the  equilibrium  limit,  in  a  similar  fashion,  by 
requiring  reaction  1,  for  example,  to  tend  toward  equilibrium  while 
reactions  2  and  3  are  unaltered.  In  other  words,  we  require  that 


^,1 


00 


(II-14) 


which  implies  that 


0^  -*  0  ,  (and  0»  -  0  )  . 


(n-15) 


Substituting  equations  (lI-6)  into  (II-9)  and  then  simplifying,  wo  obtain 


K 


1 


/  n*+  Un*+  9n»  \  /  ,  ,  \ 

(njn*.  „.n..  »•„.)  *  5|)  * 


(Il-l6a) 
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K 


2 


+ 


^  +  0(1) 


(Il-l6b) 


In  the  limit  as  tends  to  zero,  ki  tends  toward  a  finite  positive 

value  while  K2  tends  toward  +  «>  .  Equation  (ll-l6b)  in  conjunction 
with  the  rate  equation  (20)  for  q2  ,  then  Implies  that  q2  -  <l2,e  tend 
to  zero,  or 


®21<"r  "l,e>  *  <"2-  ° 


(11-17) 


To  better  understand  the  above  expression,  we  first  compute  B21 
and  via  equations  (ll-lOa)  and  (ll-ll).  When  this  is  done,  the 
following  result  is  obtained: 


®12  ■* 


®21- 


-  1 


/n*  X  /  n*  +  3n*  > 


as 


(II-18) 


Next,  substitute  equations  (l-5)  into  equation  (l3b)  thereby  resulting 
in 


^  (v  "i,e)  '  r=l,...,N2  .  (ll-19a) 

£1  ^ 

With  r  =  1  ,  with  (nj  -nj  g)  eliminated  via  equation  (ll-2b),  and  by 
means  of  expressions  (ll-loj  we  finally  obtain 

h  -  (4  *  t  "l.e’  *  '”2-  "s.e'}  ' 

in  the  limit  as  6*  tends  to  zero.  Equation  ( II -19b)  is  the  linearized 
form  of  (see  equations  (lI-3))  with  conservation  of  components 

accounted  for.  By  comparing  expressions  (ll-l?)  with  equation  (ll-19b), 
we  observe  that  the  requirement  that  0*  tend  to  zero  is  equivalent 
to  the  requirement  that  the  linearized  form  of  Lj^  tend  to  zero.  The 
equilibrium  limit  for  reaction  1  thus  reduces  the  system  of  two  rate 
equations  to  a  system  consisting  of  one  rate  equation  for  q-i  ,  discussed 
below,  and  one  algebraic  equation  fnat  represents  the  linearized  form 
of  the  law  of  mass  action  for  reaction  1. 
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When  Big  >  given  by  (II -l8),  is  substituted  into  the  transfor¬ 
mation  (II-TJ,  we  obtain  the  equation 

^1  *  -  tig  .  (11-20) 

To  better  understand  this  equation^  note  that  Li/Si  >  which  is  Inde- 
terminant  in  the  equilibrium  limit  for  reaction  1,  disappears  when  the 
difference  of  the  two  rate  equations  (II-U)  is  taken  as  follows: 


^  ^  ^  Jb. 

Dt  "Dt  Dt'  ^Sg 


(ll-21a) 


After  the  right-hand  side  of  equation  (ll-21a)  is  linearized,  this 
equation  can  be  shown  to  reduce  to 


where  Ki  is  given  by  equation  (ll-l6a)  with  6*  =  0  ,  and  equation 
(II-20)  is  used  to  obtain  q^  and  q^  g  .  Thus,  the  one  applicable 
rate  equation  (ll-21b),  in  the  equilibfium  limit  for  reaction  1,  is 
formed  by  eliminating  the  indeterminant  form  from  the  original 

rate  equations  (TI-4). 

The  foregoing  illustrates  that  the  theory  of  Chapter  3.0  has  the 
correct  behavior  in  the  frozen  and  equilibrium  limits. 
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